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I. INTRODUCTION 


А. IMPORTANCE OF DIGITAL SIGNAL PROCESSING 

One of the most pervasive advances in technology today 
has been in the area of digital signal processing. Since 
the time when the transistor made the proliferation of general 
purpose digital computers economically feasible, other ad- 
vances, particularly in integrated circuit technology, made 
practical the development of dedicated digital processors. 
From the common telephone and television, to rapid transit 
systems, manufacturing processes, speech therapy, and the 
exploration of space, digital signal processors have had 
their impact. In the Navy, digital signal processors are 
becoming more and more widely used in radar; sonar; missile 
and torpedo guidance; launcher control systems; combat 
information collection, dissemination and display; and of 


course communications systems. 


В. AREA OF INVESTIGATION 

A digital filter is defined [1] as a computational pro- 
cess or algorithm by which a digital (discrete-time and 
discrete-value) signal or sequence of numbers termed the 
МИ is transformed into an output digital signal. А 
digital filter may be implemented as a subroutine in a 
general purpose computer or as hardware in the form of a 


Special purpose digital processcr. 
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Digital filters are a subclass of the set of discrete- 
time filters which also includes "sampled data" filters in 
which the signals may take on a continuum of values. The 
mathematical treatment of the discrete-time nature of these 
filters is common to both digital and sampled data filters. 
But in digital filters the discrete value property causes 
deleterious effects in the processing of signals. 

The reason that digital filters must operate on discrete 
valued signals is that there is only a finite number of 
digits available with which to represent the signal values. 
The only way to transform one set of numbers into another 
set of numbers is by a combination of the operations of 
multiplication and addition and it is also necessary in a 
digital filter that multiplier coefficients be representable 
in a finite number of digits. In the process of performing 
a multiplication, the length of the mantissa in the represen- 
tation of the product can be as large as the sum of the 
lengths of the mantissas of the multiplier and multiplicand. 
But in general that means the product is not representable 
by the number of digits available. Therefore, some method 
of quantizing must be performed during each sample time in 
order to maintain the representability of the signal quan- 
tities. This causes non-linear effects even in a filter 
which, when viewed only as a discrete-time filter, may 
appear to be linear in nature. 

Many authors have claimed that the errors generated by 


the quantization performed in the filter are uncorrelated 


5 





ог at most negligibly correlated from sample-time to sample- 
time and between the sources of error. But when the phe- 
nomenon of limit cycles is considered it becomes obvious 
[Appendix A] that there is a great deal of correlation 
involved. Since the amount of correlation depends on the 
value of the coefficients [Appendix B] and the method of 
quantization, it follows that if there are more than one 
structure (algorithm) by which a particular filter may be 
realized and the values of the coesiicientse sand the method 
of quantization are different between these structures, then 
the errors generated and their correlation will be different. 
The objective of this work is two-foid: (1) To formulate 
a general approach for investigating noise generation and 
propagation to inelude correlation and structural erfects; 
and (2) to apply this approach to specific examples, partic- 
ularly the class of canonic second order filter structures 


which are derived herein. 


ER РПЕУТЕН ОР RESULTS 

After a summary of discrete time signal processing theory 
in Chapter II, an exhaustive search is made in Chapter III 
for all possible realizations of a second order digital 
filter in state variable form which are considered best in 
a certain sense. This search culminates in the conclusion 
that there are 58 ways to realize a given transfer function 
which meet the specified criteria. Thirty-eight of these 


realizations have not appeared previously in the literature. 
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In Chapter IV, following a method presented by Parker 
and Yakowitz [2,3], but without assuming uncorrelated errors, 
a general expression for the output error variance for any 
order filter is derived. This expression is then applied 
to the second order case and then to the realizations pre- 
sented in Chapter III. When the assumptions of uncorrelated 
and highly correlated noise are made, the expressions for 
the output error variance reduce to two sets of formulae 
for these extremal cases. Examination of these formulae 
reveals or the first time the effects of structure on noise 
generation due to quantization. Appendix D verifies that 
these formulae yield an accurate expression for the output 
error properties. 

The quantization noise problem is particularly serious 
in recursive digital filters wherein the algorithm uses the 
results of previous calculations to generate present signal 
cibles. The fact that quantization errors are also fed 
back is the root cause of effects such as limit cycle genera- 
tion. On the other hand, non-recursive filters do not exhibit 
this effect; nor do recursive filters which produce no quan- 
tization errors which are fed back. In examining the nature 
of non-recursive and recursive finite impulse response (FIR) 
discrete-time filters in Chapter V, some interesting proper- 
ties Of the signal quantities and algorithms are observed. 
These observations lead to an algorithm for computing con- 
volutions which is most efficient (in terms of number of 


operations performed) when completely overlapping blocks of 


Ши 





data from a continual cata stream need to be processed in 
this way. Additionally a new method for computing the Dis- 
crete Fourier Transform (DFT) is presented which has advan- 
tages over other methods when a continual data stream is to 
be processed. These algorithms are performed without 


recursive error generation. 
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II. DIGITAL SIGNAL PROCESSING 


A. DISCRETE-TIME PROCESSING? 


The most useful and most used tool in discrete-time 
Signal processing is the Z-Transform. Even the Discrete 
Fourier Transform is a special case of the Z-Transfornm. 
This section reviews the applications of the Z-Transforn, 
while the next section deals with the problems which arise 
when this theory is applied to digital signal processina. © 

1. The Z-Transform 

a. Definition 
Given a sequence ООН... the Z-Transform of 
the sequence is defined as 


ЕО ROT 


n- ~ СО 


» (210) 


This definition is referred to as the two-sided Z-Transform. 
Another form of the Z-Transform is defined when x(n) = 0 


for n < 0. Then (2.1) becomes 


x(n)z”" (212) 
0 


Х(2) А 
п 


Ima 


Imne use of the phrase discrete-time rather than digital 
is to emphasize that the theory in this section applies to 
"Sampled data" signal processing, in which the data can take 
Ora continuum of values, as well as digital signal process- 
ing, in which the data is restricted to discrete values [1]. 


“The object is not to reiterate known theory but rather 
to establish notation and definitions for use in later sections. 
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This definition is called the one-síded Z-Transform. Since 
(2.1) is more general, the term Z-Transform will refer to 
(2.1) unless otherwise specified. 

Both (2.1) and (2.2) can be viewed as Laurent 
series in the complex variable z. It can be shown [4] that 
for physical systems of interest (2.1) and (2.2) converge 
on the unit circle in the z-plane. Therefore, since the 


coefficients of a Laurent series are given by the integral 


m 


n dz 
27) Фе meg) Zz a (2.3) 


x(n) = 
where C is a closed contour in the region of convergence 
and encircling the origin, (2.3) is taken to be the inverse 
¿-Transform and C is taken to be the unit circle. 

b. Relation to the Laplace-Fourier Transform 

When the sequence {x(n)} is the result of the 

equally spaced sampling of a piecewise continuous signal 
then the Z-Transform, X(z) is related to the two-sided 


Laplace Transform, X(s), of the function, x(t), by 


Me )2тк 
X(z) = г Xu ea) (2.4) 
where j = М-1 and T is the spacing between samples. Note 
that this relationship is not one-to-one since there are an 
infinite number of functions x(t) which can yield the same 
set of samples x(n). Mathematically, the mapping of the 


S-plane to the z-plane by 2 = SEES many-to-one. Since 
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271 = = z+ is called the unit delay operator by the 


shifting theorem of Laplace Transform theory. 

Since the Fourier Transform can be equivalently 
viewed as the two-sided Laplace Transform evaluated at $ = jw 
then the Z-Transform evaluated on the unit circle is related 


to the Fourier Transform by 


X(z) јот ^ 
7-е 


531 
t4 8 


ХО? Е)) =4 E Xljlurku,)) 
k=-00 Паша 
(2.5) 


иПеге о is the sampling frequency in radians. Notice that 
н(е) is periodic in w with period E 
2e runet ional Transforms 

In continuous linear theory, the relationship between 


the input, u(t), and the output, v(t),is given by 


V(s) = H(s)U(s) (2.6a) 
and 
CD = J ет ат = Л а (2.6Ъ) 


со ~ СО 


where V(s) and U(s) are the Laplace Transforms of CT and 
169 respectively, ІШ is the Laplace Transform of nee 
The function h(t) is called the impulse response of the 
filter or equivalently the Green's function of the linear 


integro-differential operator which relates v(t) to u(t). 
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Equation (2.6b) is the convolution of h(t) with u(t) denoted 
by CHANG ЕС. is called the transfer function or fre- 
quency response of the filter. 

In discrete-time linear filter theory, the 


analogous equations hold, i.e. 


V(z) = H(z)U(z) (2.7a) 
v(n) = ХУ h(m)u(n-m) = X h(n-n)u(m) | (215) 
m= =% m= =% 


where V(z), H(z) and U(z) are the Z-Transforms of vín), 
h(n) and u(n) respectively. Тһе sequenceíh(n)) is called 
the impulse or unit pulse response and H(z) is the transfer 
Mmetion of the filter. 

To be physically realizable h(t) and ACD must be 
zero for t < 0 and n < O respectively. This realizability 
condition is necessary in order that the filter not react 
rozsscmetching which has not yet occurred. Thus this eendition 
is referred to as causality. 

The transfer function, H(z), can always be written 
as a ratio of polynomials in z, and H(s) as a ratio of poly- 
nomials in s. 

In designing a discrete-time filter, it is often 
the case that one tries to duplicate a well-known continuous-— 
time filter. It is then possible to transform the transfer 
function in the s-domain into a transfer function in the 
z-domain. Several functional transform methods have been 


used [5]. 
22 





a. Standard (Impulse Invariant) Z-Transform 

The Standard Z-Transform maps the transfer func- 
tion from the s-domain to the z-domain by the transformation 
2 = два the value of T and thus ш. having previously been 
decided. This method always results in the unit pulse re- 
SEANSECOI der dasemete-time filter being identical to the 
samples of the continuous-time filter taken at intervals 
T of the independent variable, t, starting at t = 0, hence 
the equivalent name, Impulse Invariant Z-Transform. There 
e ai i culty in this method; however: By Equation (2.5) 
it can be seen that the frequency response of the discrete 


КОШОК сег for -ш /2 < ш S$ wife will involve the sum of 


а 
the values of the frequency response of the continuous-time 
filter at w = © + ku» -> < К < ©, This effect is called 
aliasing and will result in H (ed 9T) not being equal to Н(јш) 
in the region of interest. They may be approximately eqval 
ШЕ ШІ! is very small outside half the sampling frequency, 
but no time limited function (which an impulse response is) 
can be limited in frequency. On the other hand, when E ju) 
is not small outside half the sampling frequency, aeu 
will bear no resemblance to the desired frequency response. 
Bor this reason other functional transforms are needed. 
bee The Bilinear Z-Transform 
Another method of transforming a transfer function 


from the s-domain to the z-domain is by the mapping 


Се (Л + 58) /(1 - 58) which has the inverse mapring 
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5 = E . This mapping, however, distorts the frequency 


response. Therefore it is necessary to counter-warp the 
desired radian frequency response befcre applying the trans- 
formation. This still does not yield an exact equivalence 
ШӨ есін the two frequency responses because the pole-zero 
lorawıons are not matched. 

c. The Matched Z-Transform 

The Matched Z-Transform is also a Bilinear 

2-Transform using the same bilinear mapping as before. But 
the way the transfer function is modified is to replace each 


pole or zero of H(s) given by 


S-s. = 0 бора) 


Z=- 7. = 2 - е i (2.8b) 


This process effectively prewarps both the radian frequency 
response and the neper frequency response. 
d. Higher Order Transforms 

Examination of the effect of the preceding func- 
tional transform techniques reveals that the algorithm implied 
by the transfer function in the z-domain is a numerical inte- 
gration routine where, for the Standard Z-Transform, the 
routine is of zero-order (rectangular) and for the Bilinear 


Z-Transform it is first-order (trapezoidal). It is possible 
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then to hypothesize using higher order numerical integration 
algorithms for the transformation between the s- and z- 
domains. This would have the effect of raising the order 
of the transfer function (which is the order of the denom- 
inator polynomial) in the z-domain. This increases the 
complexity and cost of the filter and thus is not usually 
Gone explicitly in practice, but may be done implicitly as 
observed in the next section. | 
е. Finite Impulse Response (FIR) Filters 

When a filter's impulse response (which is zero 
for n < 0) is also zero for all samples after some sample 
h(N), the filter is called a finite duration impulse response 
or simply finite impulse response (FIR) filter. In this 
case, the denominator polynomial of H(z) is zt and the 
numerator is i IE 


n=0 
to approximate a continuous-time filter of order less than 


N-n tr an FIR filter is being used 


N, then some combination of higher order integration routines 
are implicitly being performed. But the intent in designing 
an FIR filter is usually to approximate some property of 

the desired frequency response such as linear phase [6], 
minimum phase [7], or magnitude response [8,9] often without 
reference to the pole-zero locations of a continuous-time 
filter and letting the order of the filter be that which is 
needed to meet some design criterion such as minimax or 


mean-square difference. 
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3. Transfer Function Realizations 

Once the functional form of the filter has been 
decided upon, it is then necessary to design the realization 
of the algorithm to be performed. | 

a. Definitive Forms 

The most direct way to realize a transfer function 

is by performing the operations imolied by that function. 
(This is called the direct realization for obvious reasons.) 
The transfer function is usually written as the ratio of 


polynomials in z, Such as 


Ф 
N 
He 


Не еба _ (2.9) 


<, 
ним z 
e E 
о 
се, p 
N 
te 


with N > M for realizability.? Then it is possible to 


multiply numerator and denominator by 27“ to yield 


M N . 
У aa а: a 
_ 120 _ јем-М 7 _ 6:2) 
H(z) = A = ES = по) C27 О 
X E 
ј=0 7 ШЕШ 72 


37 Will be shown in Chapter III, that in some contexts 
it is necessary that N be greater than M for realizability. 
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Recalling that 271 is the unit delay operator, 


the product z^ FX(z) is the Z-Transform of x(n-k). Therefore, 
by cross multiplying the last equality in (2.10), taking 

the inverse Z-Transform of both sides term by term and 
letting Dy be unity (without loss of generality), the time 


domain equivalent of (2.10) is given by 


v(n) = 


N 
2 
1 


| са 


ay, ,u(n-1) - 


b v(n-j) (SI) 
N-M іші 5-2 


Equation (2.11) is realized directly in Figure (2-1) which 


has 2N unit delays indicated by 271, But by defining an 


intermediate variable x(n), such that 


x(n) » u(n) - 
j 


| тла 


5 Io» 


it can be shown that (2.11) becomes 


v(n) + 
1 


а 
N-M 


и м ‘2 


м3 00-1) (2.13) 
Equations (2.12) and (2.13) represent the algorithm of 
Figure (2-2) which is known as a canonic realization. It 

is canonic in that it has the least possible number of delay 


elements (N). 


b. Reduction to Lower Order Forms 
When the numerator and denominator polynomials 
ИСЕ помп in factored form, it 1s possible to realize the 
transfer function as a combination of lower order transfer 


Punctions. 
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(1) Parallel Realization by Partial Fraction 
Expansion. If the transfer function has simple poles, then 
it can be written as the sum of first and second order trans- 
fer funetions (sometimes called sections) by performing a 
partial fraction expansion. These can then be realized 
individually and then placed in parallel between the input 
and output as in Figure 2-3. If the transfer function has 
multiple poles, higher order sections will be required and 
will generally result in a higher number of delays and 
multipliers than in the canonic realization. 

(2) Cascade Realization by Factorization. In 
the factored form there is a definitive way to realize the 
overall transfer function by arbitrarily associating zeroes 
and poles and writing the transfer function as a product of 
lower (usually first and second) order sections as in Figure 
2-4, 

(3) Hybrid Realization. The operations of (1) 
and (2) above may also be used in many different combinations 
во form many different hybrid realizations, such as in 
Figures 2-5a and 2-5b. 

с. Canonic First and Second Order Sections 

Since any transfer function can be realized as 
a combination of first and second order sections, these low 
order forms take on a great deal of importance in the design 
of discrete time filters. Figures (2-6a and b) show the 
canonice realizations of first and second crder fiiters. 


Jackson [10] discussed variations of Figure (2-6b) which are 
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Just as efficient in terms of operations. Невв [11] 
extended these to 24 forms. This work extends these to 
56 forms. 
d. Other Forms of Realizations 
There are several other forms for the realization 
of a transfer functions. Among these are the so-called ladder 
structures [12,13] found by continued fraction expansions 
and the Fettweis structures [14] which are found by direct 
analogy to RLC ladder structures. Both of these types are 
ladder structures so care in terminology is necessary to 
avoid confusion. 
4. The Discrete Fourier Transform 
As mentioned in the first paragraph of this section 
the Discrete Fourier Transform (DFT) is a special case of 
the Z-Transform. The DFT is the Z-Transform of a sequence 
of length N evaluated at N equally spaced intervals around 


the unit circle in the z-plane and can thus be defined by 


N-1 : 5-1 
X(kQ) = X(z) 3kQT = E NOCERE =, КМ 
жесе n=0 n=0 
0 < k < N-1 (2.14) 
where 
2 = 2т/МТ 

a T 
ШЕ с j2r/N 
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me inverse Discrete Fourier Transform (IDET) is given by 


N-1 D 
X(kQ)W 0 < п < N-1 (2,15) 


> 

= 

5 

— 

| 
zr 
зг | 


k=0 


Notice that if n and k are outside the limits [0,N-1] then 
(2.14) and (2.15) are periodic with period N. 

The IDFT of the product of the DFT's of two sequences 
of length N is by definition reducible to the circular con- 
volution of the two sequences. The process of circularly 
convolving two sequences of length N requires ме multipii- 
cations, while the use of the DFT to accomplish the same 
operation takes 3N° + N multiplications so it is more effi- 
cient to perform the convolution directly. There were many 
early atvempts to make the DFT more efficient, including 
the Goertzel algorithm [4]. Then Cooley and Tukey [15] 
found a way to reduce the number of operations in finding 
the DFT which has become known as the Fast Fourier Transform 
(FFT) [16-19]. 

There are many variations of the algorithm now, 
including algorithms for N a composite number of mixed radix, 
algorithms which require sorting in time (decimation in time) 
as in the Cooley-Tukey algorithm and algorithms which require 
sorting in frequency now known generically as Sande-Tukey 
algorithms. Further reductions in the number of operations 
required has been accomplished by observing symmetries and 


common factors. 
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The DFT has also been widely used in the discrete 
filter design process [8, 201. 
5. Other Transform Techniques 
Other transform techniques utilized in discrete time 
signal processing include the Chirp Z-Transform [21], Fermat 
Number Transforms [22, 23], and the Walsh-Hadamard Transform 


[24,25] and the Hilbert Transform [26]. 


B. DIGITAL PROCESSING 

Wwnen the theory of Section A is applied to digital 
processing, the discrete valued Пе OL Ne Numerical 
representation required produces nonlinear effects which 
гар (пе solutions to what otherwise would be linear 
equations. Among the problem areas of concern are analog 
to digital (A/D) conversion errors, coefficient representa- 
tion (the approximation problem) and quantization errors 
involved in arithmetic operations. 

1. A/D Conversion 

The problem of A/D conversion is not one which lends 

шосе г to solution. It ТВ Prop lemeuhich must be contended 
with, the only way of reducing it being to use more signifi- 
cant digits in the digital signal representation thus making 
this purely a cost vs. error tradeoff consideration. Bennett 
[27] has presented spectral studies of A/D conversion noise, 
mis allowing the designer to represent this error as an 
injected noise at the input of the digital processor. Often 


H-souantization level of A/D converters is much larger than 
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that of the digital processor [28]. The processor may then 
include some filtering action to reduce the effect of the 
A/D conversion errors. 

2. The Approximation Problem 

When attempting to synthesize a digital filter as 
an approximation of some required characteristic such es 
impulse response or frequency response, the determination of 
filter coefficients 1s complicated by ‘the finite word length 
representation allowed. This causes error minimization 
procedures to become iterative nonlinear routines. 

The sensitivity of digital filter algorithms to 
shifts in the ideal coefficient values has been studied 
extensively. Kaiser [29] concludes that for a direct filter 
Realization, the sensitivity of pole positions increases wien 
the order of the filter. This has been end by Knowles 
and Ber 1301. Rader and Gold [31], in studying the 
e nsitivity of second order filters, conclude that two 
coupled first order sections are less sensitive than a single 
second order section. But it should be observed that two 
coupled first order sections require more multipliers than 
canonic forms (see Chapter III) and are merely a general 
state variable realization which is not canonic. Mantey 
ES studied state variable representations and concluded 
ЕЕЕ digital lters sheuld be realimed as parallel or cas- 
cade combinations of first and second order sections. This 
EN onsistent with the conclusions of Kajser above. Ladder 
МООС игез in general have been found to be less sensitive 


to coefficient inaccuracy [33]. 
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3. Quantization Errors 

Arithmetic operations in digital filters almost 
invariably lead to the requirement to requantize the result 
of the operation at some time during one cycle time. Mul- 
tiplication increases the word length of the signal. Addi- 
tion and multiplication may exceed the dynamic range of the 
numerical representation. Floating point addition may re- 
quire quantization before the operation since both the sum- 
pendsssndssddend must have the same exponent. 

duesohoriscecach ariuvchmetliesnode signilicantly etileets 
the level of quantization noise. Since fixed-point arith- 
metic is easier to implement it is most popular. But the 
relative noise level can be severe if the full dynamic range 
of the registers is not being used effectively.  Tnis requires 
sealing the driving function. Floating point arithmetic 
utilizes the exponent for automatic scaling. An intermediate 
method is known as block floating point [34] where the expon- 
ent for a whole block of data (Such as a single second order 
section in a cascade realization) is maintained in a single 
register. 

Oppenheim and Weinstein published a review [35] of 
бре effects of quantization errors in digital filter and 
Fourier Transform algorithms. This review includes the work 
oft Liu and Kaneko [36] on floating point arithmetic, Knowles 
and Edwards [37] and Gold and Rader [38] on fixed point 
arithmetic, and Weinstein and Oppenheim [39] where all three 


modes were compared on the basis of signal-to-noise ratio. 





Ер наз determined there that, in general, floating point 
Кино уде се Was best Wath block floating point falling 
between the other two modes. 

The choice of quantization mode also effects quanti- 
zation. For addition, two's complement arithmetic exhibits 
disastrous errors when overflow (upper level quantization) 
occurs [40, 41]. Zeroing and saturation arithmetic yields 
better results [42]. For multiplication and addition, 
lower level quantization methods include rounding, trunca- 
tion and sign-magnitude truncation. Jackson [43] has studied 
the interaction of upper level and lower level quantization 
errors. 

Since quantization error is a nonlinear effect, 
experience has shown that the phenomenon of limit cycies 
exists in digital systems just as in continuous ones. Limit 
cycles have been extensively studied both in statistical 
and deterministic ways. 

Certain bounds on the magnitudes of limit cycles 
in particular and accumulated roundoff errors in general in 
second order sections have been proposed by many authors and 
are tabulated in Table 2-1. 

These bounds, however, have been derived assuming 
uncorrelated error sources, a restriction which will not 


Ше assumed in the derivations of Chapter IV. 





Author Reference 
Jackson (estimate) ul 
Jackson (bound) 10 
Bonzanigo 45 
Sandberg and Kaiser 16 
Long and Trick 47 
Yakowitz and Parker 2 


f -sampling frequency 
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TABLE 2-1 Present Bounds on Quantization Error 





In statistical studies of quantization errors, the 
expressions for the output error variance, assuming uncor- 
metated errors and zero mean, uniformly distributed statistics, 
usually involve contour integrations of the form [4,48] 

2 1 


au в 6, (=) 6, (=) <Е ION 


О 2 
0 1 2T Jj 


where Ip is the output error variance 


and С. (г) is the transfer function between the zen error 


source and the output. This form of analysis requires the 
Betermination of 211 the G, (z) and the performance of all 
pne contour integrations. Furthermore, the assumption of 
uncorrelated error sources provides no insight into the 
effects of structure on the output error variance. Nor does 
this method allow the effect of assumptions made about the 


distribution and correlation of the errors to be examined. 
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Figure 2-3. Parallel Realization of H(z) 


Figure 2-4. Cascade Realization of H(z) 
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и(п) ' Ж v(n) 





H(z) = КСК, + KG, (z)6,(z) + К,6-(2)64(2)1 


Figure 2-5a. Hybrid Realization cf H(z) 





H(z) = KK, + H(z) + Н2(2) ГК, + Ha(z) + На(2) 1 


Figure 2-5b. Hybrid Realization of H(z) 


40 








Zz + a 


H(z) = ъс 


Figure 2-6a. Canonic First Order Section 


u(n) 





2 
г +а2 +a 


Н = 
(2) 2 


+ biz + by 


Figure 2-6b. Canonic Second Order Section 
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III. CANONIC REALIZATIONS OF SECOND ORDER SECTIONS 


EE INTRODUCTION 

In section IIA3b it was shown that any rational transfer 
function could be expressed as a parallel, cascade or hybrid 
parallel-cascade combination of first- and second-order 
transfer functions (or sections). For reasons of noise 
Pemeration, coefficient sensitivity, economics, logistics 
pea Others, these low order sections assume great importance 
in the field of digital filtering. 

Two ways of realizing a second order section (by the 
direct or the so-called canonic realization (see section IIA3c) 
have already been given. The "canonic" form has been and, 
in general, still is the way second order sections are 
realized in practice. However, when one considers the second 
order section in its state variable form, one realizes that 
there are an infinite number of solutions depending upon 
the manner in which the basis is represented in the state 
Есе. And, in fact, many of these are canonic in the sense 
of minimizing the number of operations required. This mini- 
ENZation also implies a minimum number of noise error sources 
due to quantization after multiplication since the number of 
multiplication operations is minimized. 

In this chapter, some inquiries are made into the relation- 
ship between the transfer funetion and the state and output 


equations of a second order section. These relationships 
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will then be exploited to find a group of state and output 
EME Ion Sets considered best in the sense of minimizing the 
number of quantization error sources. These sets are reduced 
to a number of "arrays" for compactness of notation. Half 

of the realizations represented by these arrays have not 
appeared in the literature to date nor has such an exhaustive 
search of allowable state equations been performed. The 

next chapter will delve into the ET TVs noise properties 

of the realizations derived in this chapter with a view 


toward selecting optimal structure. 


ри THE TRANSFER FUNCTION 
1. Basic Forms 
The most general form of the second order transfer 


шшпе стоп 15 given by 


0 
Н(2) = — E 005) (371) 
bo t bz" + аа Оба 


where v(n) is the output and u(n) the input. 

In this form, it appears that there are six coeffi- 
cients and hence six sources of error due to quantization 
after са оп. It is also possible to write equation 
EL) in the form | 


H(z) = (ака) 
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where 


a - 20D (3:2a) 
a, = 21/00 (3.25) 
а = a,/b, CAS) 
а = 54/90 (724) 
b = b,/b, (3.2e) 
so that there are only five non-unity coefficients. (Note 


that it is not necessary to require the condition bo A O 
since if bo were zero the denominator of (3.1) would not be 
a second order polynomial in =, hence H{z) would not be a 


second order transfer function.) So H(z) has been written 


in (3.2) in a form that has at most five quantization error 


sources. 
From the form (3.2) it is also possible to rewrite 
H(z) as 
| ae, e 
H(z) = К а * ЕЕ > В 
т ағ 2807 
where 
1 an AO ДЕ а #0 
e z (3.3a) 
0 а = 0 0 а = 0 
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a a b 
1 1 1 
— -a а #0 — - == a, #0 
Den ао 0 _ a) bo 0 
| 84 (3.3b) 
ET “0 = 0 Б. aa = 0 
0 
а а Ы 
2 2 2 
а ОИ === a, # 0 
a = MO 1 = 80 bo 0 
а, (3.3с) 
52 Еа b. E 
0 
а а #0 Е == 27 
xoc 0 0 NL 0 
(3.34) 
1 а = 0 ЛІ а = 0 


by dividing the denominator of (3.2) into the numerator in 
ascending powers of 271, 
Another form is possible if b ( or equivalently 5.) 


#0. Then the transfer function may be written 


аға! -1 
H(z) = K'a + и (3.4) 
Tot az ub 
where 
О, 
L 2 #0 Ji a, > 0 
a = = (3.4a) 
0 > = 0 0 а. = 0 


45 
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E» 
w 
CX 








EA 2 С к 
| E 
а = (ша; = 
0 ud bo зрее 
(3.4b) 
a.b a.b b 
a яо | af 
тее 2 E 20 0 
P 
"i = = 0 D а; = 0 
(з. е) 


(3.44) 


by dividing the denominator of (3.2) into the numerator in 
descending powers of 271. Note that the condition b # 0 
is equivalent to saying that there is no pole of H(z) at 
ог Елп in the z-plane. Jf there is a pole at the origin, 
then tne other pole is real and two first-order sections are 
commended rather than one second-rder Section. oO; - in 
any practical second-order section it will be considered 
ене that neither Bo пеш b is zeron 

It is still true, in the forms of (3.3) and (3.4), 


that H(z) can be realized with at most five quantization 
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EDLPOP sources. However, the coefficients K and K' need not 
be considered part of the transfer function since they are 
Peeling factors. Most often, in the practical realization 
of higher order transfer functions by cascade, parallel or 
hybrid ne (or even in realizing a second order system) a 
E 7c constant will be used which is not K or K" but rather 
Some number which, hopefully, will prevent overflow in the 
Pacers. Any gain factor necessary to restore the system to 
some overall gain factor is then included at the final out- 
put. Considering only the terms within the parentheses of 
equations (3.3) and (3.4), then there are at most four 
Amentization errors. 

If it is the case that ay = O (implying an infinite 
zero of H(z)) or а> = 0 (impiying a zero at the origin), 
Een it is possible to factor a sealing constant out of the 
rational part of (3.3) or (3.4) respectively, further reducing 
the number of error sources to at most three (four if the 
Scale factor is included). 

Or it could be that one or more of the coefficients 
is an integer in which case there can be even fewer error 
sources. But this case is not sufficiently general to 
warrant special consideration in developing algorithms 
for broad application or mass production. 

Ти the development that follows, four cases will be 
considered general enough for inclusion. These are the two 
cases of equation (3.3) with d = 0 or 1 and the cases of 


Sama tion (3.4) with d' = 0 or l. 





It will be seen that equation (3.1) or (3.2) realized 
in the "canonic" way in section IIA3c is identical to one of 
the canonic arrays to be developed, thus showing that it is 
not an additional realization. 

2. Realizability 

There is a subtle fact regarding the physical reali- 
zability of a digital filter which is not readily apparent 
in the mathematical treatment of digital signal processing. 
Because the operations involved in performing an algorithm 
cannot be accomplished instantaneously, there must be an 
ШЫ ШОП delay between the time that the input arrives and 
the time the output is available for further processing. 
If the filter is being used at very low sampling rates com- 
pared to the operation speed of the devices, this delay may 
seem negligible but it exists nonetheless. Оп the other 
МИ 11 One iS to design a filter that is usable in a broad 
field of application it is useful to consider the highest 
possible sampling rate. Then the magnitude of this delay 
is on the order of z 5, 

When the difference equation associated with equation 


(3.2) is examined, the present output, 


v(n) + agu(n) + a, u(n-1) т а „и(п-2) - av(n-1) - bv(n-2), 
(55.0) 
is seen to occur at the same time as the present input. But, 
ENNEHe previous discussion, this is impossible in a physical 
sense. Therefore we say that equation (3.2) is not, strictly 
Speaking, realizable. 
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From (3.5) it is seen that the factor which violates 
the condition of realizability is agu(n). This implies that 
(3.2) is not realizable because the coefficient, a)» in the 
numerator does not have a delay associated with it. Another 
way of saying this is that if the numerator and denominator 
of H(z) are multiplied by the lowest power of z which will 
clear H(z) of all factors, A then the power of z multiplying 
o is equal to the highest power of z іп the denominator. 
Eherefore, in order to be realizable, any transfer function, 
H(z), written in terms of positive powers of z, must have 
a numerator polynomial of lower order than the denominator 
polynomial. If H(z) does not meet this criterion, it is 
always possible to realize z" H(z) where k is an integer 
large enough to force tne transfer function into this condi- 
tion. In the case of equation (3.2), k = 1 will make the 
transfer function realizable when Ao #0. It may be argued 
that this effectively increases the order of the transfer 
function. This is of little consequence since the time 
representation of the output will be identical except for 
a shift of k sampling intervals. 

On the other hand it is not always necessary to make 
a second order section realizable in this sense. If the 
second order section is part of a cascade realization of a 
higher order transfer function, for example, it is only 
necessary to ensure that the overall transfer function is 


realizable. 
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The question of realizability will be used in the 
analysis of possible sets of state equations in the next 


section in order to reduce the number of sets to be considered. 


C. STATE EQUATIONS 

The state equations for a second order digital filter 
consist of two first order difference equations involving 
the states and input(s) of the filter. The output equation(s) 
describes a combination of the states and input(s) which will 
be seen external to the filter as the output(s). Since a 
transfer function is a relationship between one input and 
one output, only state and output equations with one input 
and one output will be considered. 

1. Basic Forms 

Define the vector, X(n), to be the vector of states 
in a second-order digital filter at time n, so that 
x4 (n) 


х(п) = | (3.6) 
xo (n) 


Now define the state and output equations of the 


filter as: 


Ax(n-1) + Bu(?) Ge 


x(n) 


| 


С =>) t sul?) (3.75) 
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where 
А 15 а 2x2 matrix 
В = 2х: matrix 
C ав а 1x2 matrix 


6 is a 1xl matrix or scalar. 


The question marks in parentheses indicate an uncertainty, 
at this point, of the proper time index to be.used. The 
n on the left side of equations (3.7) has been established 
as a reference point from which to determine the other time 
indices. Since it is already known that equation (3.7a) is 
а set of first order difference equations, it follows that 
A must multiply x(n-1) thus establishing that time index. 
Suppose the time index of u in (3.7a) is n-p, that 
EEUU (3.7b) 15 n-q and of u in (3.75), n-r, where p, а, 
and r are integers greater than or equal to zero (since 


indices implying future data are not allowed). Then 


x(n) » Ax(n-1) * Bu(n-p) (3.8a) 


Cx(n-q) * Suln-r) (3.8b) 


vín) 


Taking the Z-transform of equation (3.88) and solving 


mee X(z) yields 


Х(2) = [I - Ba Р” (3.9) 





where I is the identity matirx. Taking the Z-transform of 


equation (3.85) and substituting (3.9) for X(z) yields 


V(z) сх(2)274 * 6U(z)z | 


(C[I-Az 23 Тва (Ра) + б; Г) (г) (3.10) 


The factor inside the braces in (3.10) is thus a 
transfer function between U(z) and V(z), which can be related 
to the transfer function in question. Two things need to 
be determined, (1) the values of p, q and r and (2) the 
entries in the matrices A, B, C and ô. The first shall be 


taken up here; the second will be determined in the next 


section. 
The quantity c[I-Az i] ip in (3.10) is a scalar 
formula in z+ given by 
- a -1.-1 y + 2. 
(zZ) = C[ i-Az ] В + от Ер (Sa) 
1+ 92 + BZ 

where 

Пос b,c, ln 


e = а 20185 + ар СВ. _ 2110202 - а С. 54 (3.11b) 


a = -(а 1 + а) 





В 7 8441855 7 815854 


(Note, b, and b, here are the entries of the matrix В, not 
the denominator coefficients of equation (3.1).) If equations 
(3.8) and thus equation (3.10) are to correspond to the trans- 


her function H(z), then H(z) has to be of the form, 


H(z) 6z | 4 zz I DEC 


1) „-(р+а) 


-г (Y+ez 
87, + 7 5 


1 + ағ 


+ Bz 


a ЫН) E 852 (247) m ва) E coc Ups 
16 


ба ТГ + а 


1 + ар + bz 


(3:129 
If the requirement is now placed on the numerator 
ОГ Н(2) со be at most quadratic in nature (i.e., disregarding 
overall delay, the difference between the highest and lowest 
power of zu ШОШ О гс ег Шап two) 1t 15 possible to restrict 
the allowed values of p, q and r. 
First of all, the only indices available in the 
вере ог the filter are n and n-1, so q can only be 0 ог 1. 
Now the nighest power of zu in (3.12) is either 
ptqtl or rt2 and the lowest is either r or ptq. The following 
Loseieal arguments will place further conditions on p and r. 
a. argument (A) 
НИ ға апа рта < r then r-gtl = р ze 


EMEN Dis is a contradiction, so this case can never arise. 





b. argument B 
If ptq*tl > г%2 апа г < рға, then the stronger 
Eendrtion is that p > r-qtl, but the quadratic requirement, 
(р+а+1) - r < 2, implies that p < r-q+1, so under these 
conditions p = r-qtl. Then for q = 0, p = rtl, and for 
пи, рег, 
с. argument C 
If rr2 2 HTC DIGG] <r, then the stronger 
condition is that p ТІ but the quadratic requirement, 
wa) = (pta) < 2, implies that p ^? r-q, so that p must be 


ШКОЛЫ ӘТ р = г-да. Then for q = 0, p - г, апа fer q - 1, 


Б - г-1. 
а. argument D 
Irto Уртага ер раа слева < ри саса 
сте quadratic requirement, (rt2) - г < 2 15 automaticaliy 


КО еа. Then for q = 0, p = r or rtl, and for q= 1, 
p^ ror r-l. 

Arguments B through D taken together imply that p 
Ediffer by at most l. Recalling that p and r are 
associated with both indices of u and only u ín equations 
(3.8), it can be assumed without loss of generality that 
p and r take on only the values 0 and 1 since if they are 
greater than that, the signals u(n-p) and u(n-r) can be 
considered delayed versions of some other signal u'(n) or 
u'(n-1). Then the signal u'(n) or u'(n-1) could be used 


шансове input. 
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Teble 3-1 lists all possible combinations of p, q 
сит сисп that they take on only the values 0 and 1, along 
with the associated transfer function. Іп addition,the 
ШӘРПев ог ptq, ptqtl and rt2 are tabulated along with the 
appropriate logical argument which applies to their relative 
sizes. Referring to those arguments and the corresponding 
values of p, q and r, it is seen that lines (b) and (g) 
of Table 3-1 violate the conclusions of the argument with 
respect to the relative values of p and r. Thus the trans- 
fer functions of lines (b) and (g) cannot be related to a 
second order transfer function whose numerator is at most 
quadratic in nature. 

The state equations and output equation of line(h) 
ши арте 3-1 can be seen to be identical to those of line 
(c) except that the input in both equations of (h) has been 
delayed one sample period. Therefore the internal structure 
would be the same except for a delay in the input line. 

The equations of line (e) can also be seen to be identical 
in nature to line (h) by rewriting the output equation of (e) 
as v(n-1) = Cx(n-1) + óu(n-1) instead of 

v(n) = Cx(n) + Su(n), whereas in (h), v(n) = Cx(n-1) + Suln-1). 
in other words, the output of the equations for line (h) 

is a delayed version of the output corresponding to line (e), 
So the structure will be the same except for a delay in the 
Gm@cDut iine. 

Similarly the equations of (f) represent a delay 
added on the input line of (a) and the equations of (d) 


represent a delay on the output line of (a). 
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These two DS of essenvially identical state and 
output equations are indicated as groups S and T respectively 
in Table 3-1. (Note that for ô = 0 any of the transfer 
functions in Table 3-1 could belong to either group, including 
the two that apparently violated the logical arguments since 
for 6 = 0 the value of r is irrelevant. In order to remain 
as general as possible, ô will remain arbitrary.) Thus the 
transfer functions associated with each member within a 
particular group represent essentially the same transfer 
function with the possible exception of an overall delay. 

For simplicity, two transfer functions are chosen to 


represent their respective group, namely 


на(а) = [6 + Фб(аја ја (3.13) 


апа 


H(z) = 161 + (2) 1271 cam 
Where ©'(z) has the same form as %(2). Equation (3.13) is 
taken from line (h) of Table 3-1 and (3.14) from line (f). 
The primes in (3.14) have been added for purposes which will 
become clear in the next section. The reason for the group 
designations and subscripts, S and T, is historical and will 
be made clear at the end of this section. These two forms 
have peen chosen because dias: an overall delay and are 


thus physically realizable. 


(| 





It is now possible to make some specific connections 
between the two forms of H(z) given in equations (3.3) and 
(3.4) and the state equations associated with (3.13) and 


(3.14) respectively. 


р. THE COEFFICIENTS OF THE STATE EQUATIONS 

Given a transfer function in the form of equation (3.3) 
(with an overall delay added for realizability) and a set of 
state and output equations which yield a transfer function 
in the form of (3.13), the two systems can be made equal by 
setting their transfer functions equal and solving for the 
appropriate coefficients.! In this case (disregarding 


Scaling constants) 


| 
е 
E 
Q 
N 


| -1 -2 \ 
H(z) -1 р ET ee, ен 


із az + 5277] 


Не(а) (3.15) 


ExuSeing coefficients, it is found that for (3.15) to be 


compatible 


cuc = -(a44 + 270) (3.164) 


тре derivation which foliows is essentially due to Parker 
and Hess [49] which follows a similar effort by Jackson L10,. 
Meenderivation is reiterated here for continuity. The only 
difference is some minor changes in notation and the use of 
 атар1е forms of the transfer function. Parker and Hess 
used equation (3.3) without an cverall delay and the State 
equations corresponding to line (e) of Table 3-1 with & change 
NN all time indices of the state equations but no change 


E e output equation. 
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о = В 814855 7 845851 (3.16b) 
с = у = сур, + cob. (3.16c) 
а = 6 (20059 


2126162 * 8540€5b, - 844106555 - 8550154 


(3.16e) 
Ерпс that the coeificients а, b, c, d and e are given, 
it is necessary to solve equations (3.16) for the entries 
of the matrices of the state and output equations from which 
(3.13) was derived. Any combination of the ajj? b. and c, 
which yields the given values of a, b, c and e is a valid 
realization. But there are eight coefficients in the matrices 
AS B and C which implies that in some realizations there 
could be as many as eight quantization error sources. Further-- 
more if the values of a, b, c and e are given to k-bit accur- 
acy, the solutions to equations (3.16) will not, in general, 
yield coefficients to k-bit accuracy. The approach taken 
by Parker and Hess was to force the solutions to be in k-bit 
accuracy by writing (3.16a, b, c and e) as 


== (аз жал) (3.17а) 


11 22 


а 1222 4 - 1942221 4 (3.17b) 


2d 





© 
! 


(516114 Е [b,c a (3.170) 


ER | 


q b, 


47 181199244 7 [2559454 Jq 


(3.184) 


where [ * ] d indicates quantization of the product within. 
Since there are four equations in eight unknowns, it is 
possible to choose an appropriate set of four unknowns arbi- 
marily and solve for the other four. Furthermore, the 
solution is simplified by the intuitive realization that 
quantization of product terms can be avoided if, in the 
products of (3.17), each term except one is unity or one of 
eem is zero. Coincidentally, it turns out that this approach 
Beso results in only fleur non-zero, non-unity coefficients, 
1.е., the four coefficients remaining after arbitrarily 
Eoszunge four to be unity or zero. This approach thus 
yields the minimum number of quantization error sources as 
pell. 

By applying a set of rules for the selection of 
unity and zero coefficients, Parker and Hess obtained eight 
system март сев“ which would meet the criteria. Only two 


of these were unique, two others being transposes of the 


А: В 

Exnsse system matrices were of the rom [LL 628 
ИСО опїпе the matrices yields the С. 

Шо сез А, B, C and 9. The phrase "system 

array'"will be preferred here for reasons to be shown later. 
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first two, the other four simply involved the interchange 
of rows or columns in the first four. 
The two unique system arrays derived by Parker and 


Hess were 


-а -b 1 
= о (3.18а) 
с е 
-(1+а) -(1+a+b) I 
S, = 1 1 0 (3.100 
e cte d 


The 5, array had been previously discussed by Jackson. The 


transposes of these arrays are then 


r-a T с 1 
= in 0 e | (3.192) 
1 0 а 
| -(1+a) 1 с 
5 - | -(1tatb) 1 cte (3.195) 
1 о а 


Jackson showed that the transpose of a system array 
yields a configuration identical to the untransposed config- 
Merion except that (1) the direction of signal flow is 
reversed and (2) summing nodes become pickoff nodes and vice 
versa. Jackson further showed that for each unique system 
бау апа its transpose, twelve configurations can be achieved 
through flow graph manipulations. The conclusion of Parker 
and Hess was then that there were twenty-four canonic 


Pealizations. 
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The previous analysis was based on equating the 
transfer function of equation (3.3) with the transfer function 
of (3.13). A similar analysis can be carried out by equating 
(3.4) to (3.14). Then 


EX -1 
! - 
H(z) = (а! + Ба! 1 = (6! + Е 
jut Ча t Diz | l + ag + BZ 
= На( 2) | (3: 200 


шише а", е!, 6!, Үү! апа є' аге ргіпей to distinguish 
them from the analysis which precedes, while a, b, œa and $ 
are not primed since the characteristic polynomial must be 
the same for all realizations. When the coefficients are 
Pemaved, however, it is seen that the equations are identical 
in form to equations (3.16). Thus following the same set 
шине for selection of unity and zero coefficients, the 
system arrays which result will have the identical form of 


5. and S, of equations (3.18). The only exception is that 


b 
Will be replaced by d', e by c', and e by e'!, which means 
that for the same transfer function H(z), the arrays which 
result from (3.20) will differ from those resulting from 
(3.19) numerically. In order to indicate which form of the 
Geemoter function is used to find the entries in the array, 
а Т will be used for those resulting from equation (3.4) 
whereas an S designates the ones found by Parker and Hess, 


mao the rationale for calling the groups 5 and T comes 


PmeomenisvOorical development. The next section wili force 
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another change in notation (altering the subscripts, e.g. ‘a’ 
in 5.) to avoid an historical conflict and the resulting 
confusion that might be created. 

In summary, the state equations and associated 


transfer functions are given by 


х(п) к Ах(п-1) * Bu(n-1) 
>Н(2) - На (2) 
v(n) = Cx(n-1) + би(1-1) 


Sa cz ^ * eg"? 2 
B ST Ten 
1 Реза 42 1017 
(3.21) 
апа 
x(n) = Ax(n-1) + B'u(n-1) 
y *H(z) = Hp(z) 
v(n) = C'x(n) + $'u(n-1) 
-1 
EN Мб EE s 271 
Вах ре 
(3.22) 


E. CANONIC ARRAYS 

When Parker and Hess wrote their state and output equations, 
it was done in such a way as to make it possible to write 
all three equations as one matrix equation. Thus 


mht je eax) + Buln) (3.232) 


Ша) Cn) t óuln) (3:220) 
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because 


i 
E. E Cau 21 
x(n+1 І х(п) 
" m ejje Sod?» | | = | (3.24) 
vin TE u(n 
Cy с. 6 


Thus, the arrays S, and S, of equations (ето теошја truly 
пенса са matrices. Examination of Ehe entries in Table 3-1 
shows that by raising the time indices of the state equations 
f line (e) by 1 to coincide with (3.23a), the same result 
is achieved. ТЕ 15 also possible to write the state and 
output equations of lines (c) or (h) directly as one matrix 
equation. These three lines are all in group S. When one 
tries to do the same thing with che equations in group T 
EOS round to be impossible. Therefore, the arrays asso- 
ciated with group T cannot be called matrices. Hence, in 
this work only the description "array" will be used when 
referring to the 3x3 arrays S or T. The transpose of an array 
is defined as though it were a matrix. A, B, C and 6 continue 
to play the role of matrices and will be referred to as such. 
As mentioned earlier, Jackson showed that there were 
six realizations possible for an array and six for its trans- 
pose. This can be accomplished by flow diagram manipulations 
or, equivalently, by manipulations of the equations involved. 
UNO of the six are identical except that in one d = 0 and 


in the other, d = 1. In the remaining four it was claimed 
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that d had to be unity. In fact, it will be shown that there 
are nine manipulations of the equations possible for the 
State equations of group 5 апа that each allows d to be 
zero or unity. There is a realization for all 36 of these 
cases (the nine manipulations of the array S,> nine of S., 
each with d = 0 and d = 1). Figures 3-1 through 3-18 show 
these realizations. In addition, there are the 36 transpose 
configurations of these. To allow the realization of any 
transfer function with a single realization, composite 
eanl zations are included which usually require an additional 
summer and which change the configuration from 9 = 0 tod= 1 
by the use of a single switch.? 
Furthermore, the state equations in group T can also 
be manipulated to yield different flow diagrams. There is 
not as much flexibility in this group, however, and only forty- 
four unique configurations have been found (compared to 
72 for group S). Figures 3-19 through 3-22 show the 
untransposed realizations for group T. 
Since previous work has referred to the realizations 
іп a manner which is not adaptable to a one-to-one corres- 


pondence with the realizations derived here, the notation 


needs to be changed to avoid confusion. 


~ 


A switch will be indicated as follows (see Figure 3-23): 
a) A switch which is normally closed for d = 1 and open 
for d = 0 is represented as a multiplication by d. 

b) A switch which is normally open for d = 1 and closed 
for d = 0 is shown as & multiplication by (1-4). 


3 
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Call an array in the form of equation (3.182) an 
"M" appay; an array in the form of equation (3.18b) is an 
"N" array. When the array is realized for the form of the 
transfer function and state equations of (3.21), precede 
its designation by an "S", In this case the entries of the 
array include the unprimed quantities: c, d and e. Table 
3-2a lists all the SM array manipulations. Table 3-2b lists 
those of the SN array. SM compares to the За array of 


00 


(3.18a). SNoo compares to S, 


When the array is realized for the form of the trans- 


in (3.186).4 


fer function and state equations of (3.22), precede its 
designation by a "T", In this case the entries of the array 
include the primed quantities: c', d' and e'. Table 3-3a 


lists al) the TM array manipulations. Table 3-3bd snows спе 


single array, TN This compares to the S, array with 


0” b 
primed quantities replacing the unprimed. Similarly ІМ, 
compares to Sa with primed entries. 

Мен сгоар © it will be shown that the transpose 
configurations are not the transposes of Tables (3-3) but 


rather correspond to those of Tables 3-2 by substituting 


primed quantities where applicable. 


‘the subscripting notation will be explained after an 
meamination of the process of manipulating the state 
Fouations. 
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| The development of these realizations follows. (Only 
operations which do not increase the number of non-zero, 
non-unity multiplications are allowed.) 
Figure 3-la shows the flow diagram of the array, 
5М with a = 1. The state and output equations for this 


00 
realization are, according to equation (3.21): 


x, (n) = -ax, (n-1) - bx, (n-1) + u(n-1) (3.25a) 
xj(n) = x, (n-1) ( 3.2209 
v(n) = ex, (n-1) + eX, (n-1) + лошата (0272505 


Figure 3-1b shows the realization of these equations 
Bor d - 0. This configuration corresponds to the realization, 


S al of previous work. Figure 3-1а corresponds to 5. 


2. 
Figure 3-1с exhibits the composite realization which 

allows 9 to be either 0 or 1. When Parker and Hess said there 

were 24 canonic realizations of a transfer function, the 

meme Was, that they found four realizations for a transfer 

function with d = 0 and twenty for a transfer function with 

ELE но превзе аге never the same transfer function. There 

is a class of functions with @ = 0 and a mutually exclusive 

class with d = 1. Correspondingly, there is a class of func- 

tions with d' = О апа а mutually exclusive class with d' = 1. 


Taken together, there can be overlap between the classes of 


EEUU —Ihat is if d 9» O, d' could be 0 orl, etc. In 
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total there are four classes given by (d = 0, d' = 0), 
(а= 0, а! = 1), (а = 1, а! = 0) апа (а = 1, а! = 1) which 
include all second order transfer functions. The inclusion 
of a composite realization is in the interest of logistics. 
Each composite realization can be used for any second order 
ШватвГег function. 

Proceeding with the development of the other forms, 
consider equation (3.25c). This equation can be rewritten 
by adding and subtracting the quantity bx, (n-1) to yield a 


new output equation 


v(n) = cx, (n-1) + (bte)x,(n-1) * du(n-1) - bx, (n-1) 
(5225) 
The quantity -bx,(n-1) is one which is formed in calcuiating 
x4, (n) Bm cis" avaidable to be used in (3.20 without adding 
another multiplier to the realization. 
Equation (3.26) together with (3.25a and b) are 
realized in Figure 3-2a for d = 1. In Figure 3-2b, d = 0 
ша 3-2c and d show the composite realization. These are 


realizations of SMo1 where the array is 


el 





On the other hand it is possible to rewrite equation 
(3.25a) by adding and subtracting exj(n-1) to the right hand 


side. Then 


xy (n-1) = -ax, (n-1) - (bte) x, (n-1) + une) + ex, (n-1) 
(See) 
Equation 3.27 together with equations (3.25b and c), 
the original second state equation and the original output 
equation, are realized in Figures 3-3a, b, c and d with 
а = 1, а = 0 апа two composites, respectively. In this case 


the array is SMo> given by 


-2 e-(bte) al 
1 0 0 
© е а 


It is now possible to describe the subscript notation used. 
In SMoo there are no manipulations of the original array. 
In SM and SM the second subscript is changed to indicate 


ol 02 
a change in the second colum of the array. When the sub- 
Spt is l it indicates adding and subtracting the coeffi- 
Ж ЕПС ОТ the first row in that column to the coefficient 
Ot the last row in that column. A subscript 2 indicates 
the opposite operation. These operations can be performed 


ШІ (ле first column as well, or in both columns at the same 


tame and in the same or opposite directions. A complete set 


Да 





of arrays with all these operations гге given in Tables 3-2a 
and b. Only the untransposed arrays for group S are given 

in the tables. Figures 3-1 through 3-18 show all the untrans- 
posed configurations implied by Table 3-2, whereas Figures 
3-24 through 3-41 show the realizations found from the 
transposes of the arrays of Table 3-2. 

Шера 15 the form historically called the 
Beanonic'’ realization if ag is factored out of the numerator 
of equation (3.2). 

For the arrays of group TM, manipulations are only 
possible between the two entries corresponding to 871 and 
со, that is the first entry in column 1 of the array and the 
last entry in column 2. (All manipulations in group S were 
accomplished within the same column.) There are apparentiy 
no ways to rearrange the equations of array IN, without 
шлогодистпре more multipliers and summers. 


To change TM, into TM consider the state and output 


0 JU 


equations for TM Нәсіп еацавтоп (3.22), 


gr 
x4 (n) = -ax, (n-1) - bx, (n-1) + u(n-1) (3.28a) 
EE x (n-1) (3.285) 
Ка (ап) + е'х,(п) + du(n-1) (3.280) 


By adding and subtracting the quantity ax, (n) on the 


right side of (3.23c) the output equation becomes 





v(n) = c'x (n) + (ate! )x5(n) - ах. (п) * du(n-1) (3.29) 
But by equation (3.28b), X, (n) - X, (n-1) зо 
v(n) + c'x, (n) + (ate')x,(n) - ax, (n-1) * du(n-1) (3.30) 


The quantity -ax, (n-1) is one which is formed in 
calculating x, (n) so that it is available to be used in (3.30) 
without adding another multiplier to the realization. 

Equation (3.30) along with (3.28a and b) are realized 


in Figure 3-20 ав ТМ The opposite manipulation indicated 


1: 


by TM, in Table 3-3 is shown in Figure 3-21. 


2 


The untransposed realization of TN, is given in 


Figure 3-22. 


mie ana IM. have no transpose СОГ ОО but 


1 2 


TM, foes. It is not true, however, that the transpose of TMo 


is realized by reversing the flow in Figure 3-19 and inter- 
changing summing nodes and branching nodes. The realization 
is found by applying the state equations of equation (3.22). 


Furthermore, taking sm; ери ет о екп ару ско л 


d' respectively and using the state equations of (3.22) yields 
another set of realizations for H(z). These realizations are 


designated TMs y (with i and j corresponding to the subscripts 
E ) Note that TH and тм? 
CM 0 00 


The same correspondence is true for the transpose of 


of SM are the same realization. 


Т. Taking D. псп ас1по о, е, апаа by c', e' and d* 


211 the TMj and TM2 arrays are transposed it is seen that 
Eve multipliers would be required, thus they would no longer 
ре сапопіс in terms of multipliers and error Sources. 
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respectively and using the state equations of 3.22 yields 
the realization TN, ys with TN¿, the same as TNG. 
Figures 2-42 through 3-59 show all the transpose 


configurations for the T forms. 


Е. CONCLUSIONS 

This development has included a very exhaustive search 
of allowable state equations and transfer functions for 
realizations of second order digital filters, canonic in 
the number of delay elements and multipliers (and thus quan- 
tization error sources after multiplication) and also canonic 
in that the multiplier coefficients are linear combinations 
of the coefficients once they are given in the S or T form 
Meche transfer function. 

in brief, the arrays of Tables 3-2 and 3-3 when combined 
with the state equations of (3.21) and (3.22), respectively, 
yield the realizations of Figures (3-1) through (3-22). The 
transposes of the arrays of Tables 3-2 when combined with 
both sets of state equations yield the realizations of Figures 
(3-24) through (3-41) for Equations (3.21) and Figures (3-42) 
rough (3-59) for Equations (3.22). 

In total, these realizations consist of 36 realizations 
for transfer functions with d = 1 and 36 for d = 0, as well 
as 22 realizations for d' = 1 and 22 for d' = 0. That is, 
for a particular transfer function, there are 58 realizations 
Eimce d cannot be both O and 1 nor can d' be both 0 апа 1. 


Included are composite realizations which allow d or d! to 


O 





be either 0 or 1. These require the inclusion of a switch 
to indicate the particular value and sometimes an additional 
summing device. 

The concept of "canonic arrays" introduced here made 
this development possible. If one forces the array to play 
the role of a matrix then only the state equations of the 
S forms are possible. This forces the transfer function 
also to be in the S form. The T ‘forms only arise when one 
allows for the concept of canonic array or, equivalently, 
separating the operations of the state equations from those 
of the output equations. This allows the time index of tne 
states on the right hand side of the equations to be different. 


This factor is essential in developing the T forms. 
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Figure 3-la. Realization of SMog with d=] 





Figure 3-1Ь. Realization of SMoo with d=0 
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Figure 3-1c. Composite Realization of SMag with d=0 or 1 
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Figure 3-2b. Realization of SMa] with d=0 
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Figure 3-2c. Composite Realization of Sro] with d=0 or Y 





Figure 3-2d. Composite Realization of SM with d=0 or 1 
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Figure 3-3a. Realization of Shiga with d=] 
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Figure 3-3b. Realization of SMgo with d=0 
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Figure 3-3c. Composite Realization o 
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Figure 3-3d. Composite Realization of Мо? Wich а = Оно | 
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Figure 3-4b. Realization of SM0 with d=0 
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Figure 3-5a. Realization of SM), with d=] 
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Figure 3-5b. Realization of SM,, with d=0 
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Figure 3-5c. Composite Realization of SM,, with d-0 or 1 
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Figure 3-5d. Composite Realization cf SM, y wtih 4-0 оғ 1 
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Figure 3-6a. Realization of 5М 2 with d-1 
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Figure 3-6b. Realization of Е with d=0 
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Figure 3-6c. Composite Realization of SM 2 with d=0 ог 1 
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Figure 3-6d. Composite Realization of УМ. wtih d-0 or ] 
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Figure 3-/a. Realization of SM,, wtih d=1 
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Figure 3-7b. Realization of 5,7 with d=0 
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Figure 2-7d. Composite Realization of SM o with d-0 or 1 
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Figure 3-8a. Realization of SMa, with d=] 
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Figure 3-8b. Realization of SM,, with 4-0 
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Figure 3-8d. Composite Realization of SM,, with d=0 or 1 
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Figure 3-9c. Composite Realization of УМ, with d=0 or 1 
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Figure 3-9d. Composite Realization of SH,» with d=0 or 1 
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Figure 3-10a. Realization of SNoo with d=] 


u(n-1) v(n) 


-(1+а+5) ) с+е 


Figure 3-10b. Realization of SNog with d=0 
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Figure 3-10c. Composite realization of SN.. with d=O or 1 
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Figure 3-lla. Realization of э] with d=] 
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Figure 3-1lb. Realization of SN,. with d=0 
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Figure 3-11d. Composite Realization of Shay wtih d=0 or 1 
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Figure 3-12a. Realization of SN). with d=] 
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Figure 3-12b. Realization of SNo? with d=0 
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Figure 3-12c. Composite Realization of SNo, with d=0 or | 





Figure 2-12d. Composite Realization of 516% with d=0 or 1 
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Figure 3-13a. Realization of УМ. с with d=] 
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Figure 3-13). Realization of УМ С with d=0 
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Figure 3-13c. Composite Realization of 5840 with 4-0 оғ 1 
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Figure 3-13d. Composite Realization of SN,¿ with 4-0 оғ 1 
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Figure 3-14b. Realization of S1 witn d=0 
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Figure 3-14d. Composite Realization of 5811 with d=0 or 1 
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Figure 3-16a. Realization of SNon with d=] 
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Figure 3-16b. Realization of Shing with d=0 
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Figure 3-16c. Composite Realization of SN. with d-0 or ] 





Figure 3-160. Composite Realization of SN,» with 4-0 оғ 1 
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Figure 3-17a. Realization of SNoy with d=] 
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Figure 3-17b. Realization of SN» with d=0 
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Figure 3-17c. Composite Realization of SN,, with d=0 or 1 
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Figure 3-17d. Composite Realization of SN»: with d-0 or 1 
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Figure 3-18a. Realization of SN5, with 4-1 
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Figure 3-18b. Realization of Noo with а-0 
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Figure 3-18c. Composite Realization of SN,, with (6-0 оғ 1 
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Figure 3-18d. Composite Realization of Но with d=0 or 1 
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Figure 3-19a. 
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Figure 3-19c. Composite Realization of ТМ, wtih d'=0 or 1 
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Figure 3-20a. Realization of TM, with d'=] 
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Figure 3-20b. Realization of ТМ, with d'=0 
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Figure 3-20c. Composite Realization of ТМ, with d'=0 or 1 





Figure 3-20d. Composite Realization of TM with d'=0 or 1 
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Figure 3-2la. Realization of TM, with d'=] 
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Figure 3-21b. Realization of TM, with d'=0 
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Figure 3-21c. 
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Figure 3-22a. Realization of TN, with d'=] 


v(n) 


u(n-1) E | | © 


| / 
'9 
с! 
де 
| 


-(1*a) 





-(1+atb) 


Figure 3-22b. Realization of ТА with d'=0 
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Figure 3-22c. Composite Realization of TN) with d'=0 or 1 


120 





-( « » multiplier 


— 
Switch 
AND gate 

МЕС for d=] 

GND for d=0 


(a) POSITIVE LOGIC 


multiplier 


№.0.(4=1) 
--7- switch 


puo AND gate 


GND for d=] 
ү for d=0 
GC 


(b) NEGATIVE LOGIC 


normally open 
normally closed 


a 
Сэ © 


Figure 3-23. Equivalent symbols for digital switches 
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Figure 3-24, Realization of SM, with d=0 or 1 
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Figure 3-25. Realization of SMyy with d=0 or 1 
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Figure 3-26. Realization of SM with d=0 or ] 


u(n-1) i vin) 
p — _ AE 
d=0 


x(n) 








Fiugre 3-27. Realization of SM, o with 6-0 оғ 1 
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Figure 3-28. Realization of SM, 
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Figure 3-29. Realization of S with d-0 or 1 
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Figure 3-30. Realization of 5м, С with d=0 or 1 
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Figure 3-31. Realization of SM, y with d=0 or 1 
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Figure 3-32. Realization of SM 
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Figure 3-33. Realization of Мо with d=0 or 1 
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Figure 3-35. Realization of SNo% with d=0 or 1 
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Figure 3-36. Realization of Silao with d=0 or 1 
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Figure 3-37. Realization of SN] with d=0 or 1 
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Figure 3-38. Realization of SN. 5 with d=0 or 1 
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Figure 3-41. Realization of S with d-0 or ] 
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Figure 3-43. Realization of ТИ, with d'=0 or 1 
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Figure 3-44, Realization of TM. with d'‘=0 or 1 





| v(n) 
> E > 
хо(п) 





+ 
2 


Figure 3-45, Realization of ти. ( with d'=0 ог | 
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Figure 3-46. Realization of ТМ. 
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Figure 3-47. Realization of ТМ. 5 with d'=0 or 1 
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Figure 3-49. Realization of тм): with d'=0 or 1 
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Figure 3-50. Realization of TM, with d'=0 or 1 
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Figure 3-54. Realization of TN,, with d'=0 ог 1 
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Figure 3-55. Realization of IN. with d'=0 or ] 
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Figure 3-56. Realization of TA with d'=0 or 1 
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Figure 3-57. Realization of тк, with d=0 or | 





Figure 3-58. Realization of ТМ with d'=0 or 1 
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ive TERROR ANALYSIS OF RECURSIVE DIGITAL FILTERS 


A. INTRODUCTION 

The effect of finite precision arithmetic on the perfor- 
mance of digital filters has been considered from the deter- 
ministic (limit cycle) and probabilistic (random noise) point 
Of view by several authors [3,!5-48]. Calculation of output 
noise due to quantization, such as round-off after multipli- 
moron Of Signals by constant factors, usually involves a 
difficult contour integration [48]which does not readily lend 
сет to physical interpretation in terms of filter structure 
ШЕ сопГігрпгасіоп for purposes of comparison. In addition, 
the usual noise formulas are based upon the assumption that 
the quantization error sources within the filter are uncor- 
related. Although this assumption has been confirmed as 
reasonable by several authors [35,50], it is interesting to nocte 
that under the conditions of a limit cycle, quantization 
ШЕЕОгв within the filter are highly correlated. Thus the 
lufluence of correlation on noise calculations remains a 
EDestion. 

In this chapter, several results are presented. FI SER 
following a procedure introduced by Parker and Yakowitz [3], 
but without assuming uncorrelated error sources, the covariance 
matrix of the state errors and the variance of the output 
error are expressed in terms of the correlation matrices of 


the quantization error sous This enables the effects oz 


lThis approach was presented at the Arden House Workshop 
on Digital Signal Frocessing [57]. 
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correlation between Sources to be taken into account. Тһе 
results are presented in a general form which enables the 
structure of the filter to be introduced independently. 
More specifically, the expression for the output error 
variance is applied to second order sections and then to the 
Benonie realizations of Chapter III. 

Experimental verification of the results are presented 


in Appendix D 


B. STATE SPACE FORMULATION 
In Chapter III the equations for a digital filter in 


state variable form are given by 


x(n) » Ax(n-1) * Bu(n-1) (4.1а) 


v(n) = Cx(n-1) + 6u(n-1) (4.1b) 


2 
Сог the S forms (see Chapter III, Equations (3.21)) , and 


! 


х(п) Ax(n-1) * B'u(n-1) (4.2a) 


v(n) СОЕ ó'u(n-1) (17259 





“The use of the terminology S forms (or T forms) is not 
restricted to second order filters. The state and output 
equations in matrix form are in general applicable to nth - 
Mer filters. Even to the transfer function of an nt) order 
Ша Бет, the S,T terminology applies in general. 5 means 
dividing the numerator by the denominator in ascending powers 


of 2-1. T means to divide in descending order. 
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дот сие T forms (from Equations (3.22)), where А, B, B', C, 
С', 8 апа 6' are coefficient matrices; x(n), the vector of 
State variable at time nT; u(n), the input; and vin), the 
риорцс. 

If one assumes finite precision arithmetic in carrying 


out the digital filter algorithm, then Equations (5.1) become 


ки) LAx*(n-1)1, + LBu(n-1) Jc (4. За) 


um [Cx*(n-1) J, + (646-1710 (4.35) 
where the asterisk denotes rounded signal variables and 

[ нЕ denotes the process of quantization after multiplication 
of a previously quantized signal variable by & multiplier. 


x 
Bhav is, if x; (n) is the jeh component of x and as is 


the ij° entry in Шешпей “һе КЕН component of LAx* (n-1) J. 
p + 

13 given by 2 L8, 4X, (n-1)J, where p is the order of the 
= | 

Euer. 


Following Yakowitz and Parker [ 2 |, the following bounded 


vectors can be defined: 


e(n) » Ax*(n-1) - [Ax* (n-1) Jg + Buln-1) - [Bu(n-1) J, 
(4 4a) 





SAlgorithms which quantize after summation are also 
possible aná have been studied [ |. Further discussion of 
this and other methods are considered in Section D. 
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е(п) = Сх%(п-1) - LOx*(n-1)], t $u(n-1) - [éu(n-1) ], 
(4.45) 


Ehese are the quantization errors introduced into the filter 
at each time, nT. There can be several contributions to each 
component of e(n) depending upon the number of multipliers 
in the appropriate row of A and B which require the product 
to be quantized; and similarly for e(n). 

If one defines the error, y(n), as the difference between 
the states of the finite precision realization and those of 


ЕМЕ ideal infinite precision realization, then 


ШІП) - хп) - x*(n) 


| 


Ax(n-1) * Bu(n-1) - [АХ“(п-1) Је _ [Bu(n-1) Jo 
(4,5) 


by Equations (Y.la) ana (4.3a). The error, y(n), is by 
definition the accumulated error in the states of the finite 
Precision realization. In order to find an equation Гог the 


propagation of this error, Equation (4.5) can be rewritten as 


|! 


ШОУ = Ах(п-1) + Ви(п-1) - (Ax (n-1) 7, - [Buln-1)1, 


+ Ax*(n-1) - Ax*(n-1) 


A[x(n-1) - x*(n-1)] + Ax*(n-1) - [Ax*(n-1)1, 


* Bu(n-1) - [Bu(n-1) ], 
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y(n) = Ay(n-1) + e(n) | (4. ба) 


Thus, the accumulated error in the states at time nT depends 
Encthe previously accumulated error (at time (n-1)T) and the 
Epor introduced by quantization at time nT. 

Applying Equation (4.6a) recursively and assuming that 
there is no accumulated error at time.nT - 0, i.e. x*(0) = x(0), 
еп the dependence of y(n) on all previously introduced 


quantization errors can be found to be 


n Е 
y(n) = 5 А!7%е(3) (4.6b) 
j=1 


ШОУ a change of variable in the summation, 


n=l 


ШІП) = 2 


1=0 


діе(п-і) (8.66) 


Equations (4.6a-c) show that the effect of quantization 
errors introduced "long before" time nT becomes small since 
Ши на 1 сег 15 assumed to be stable and the characteristic 
equation of the filter depends solely on the matrix A. What 
may be considered "long before" depends strongly on the 
characteristic of the fisher ND CHEERS for second 
Order filters, a low damping СЕ Тепе would imply a con- 
Siderable time for the effect of a quantization error to 
become negligible. This idea will be useful in evaluating 


the results of the derivation and will be seen to be applicable 
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бо the stochastic formulation in Section C as well as in 
the deterministic Equation (4.6). 

Defining Av(n) as the difference between the outputs or 
the finite and infinite precision realization, a similar 


derivation yields 


Ау(п) v(n) - v*(n) 


| 


Cx(n-1) + Su(n-1) - [Cx*(n-1)] - [$u(n-1)] 


+ Cx*(n-1) - Cx*(n-1) 


Cy(n-1) + e(n) 


п-1 1-3 
= C БОА 79 е(3) си ба) 
UE 
n-2 . 
= таре) (4.64) 
1150 


The quantity, Av(n) is the accumulated error at the out- 
put of the digital filter at time nT. Equation (4.64) shows 
that this accumulated error depends on all previous quanti- 
zation errors generated in the state equations (not including 
mee errors generated at time nT) but Av(n) only depends on 
the present quantization errors generated in the output 


Eguation. 
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Equations (4.6) apply to the S forms of the state eaua- 
tions. In order to do the same for the T forms it is neces- 
sary to redefine the quantization errors in the output 


equation at each time nT as 


са с Охи) - Loiz J, + S'u(n-1) - [Stu(n-1)], 
(4.7) 


Keeping the definition of e(n) the same as (4.4a) except 
to replace B by B' the accumulated state and output error 


equations become 


D n-j пеш 
y(n) » Ay(n-1) + е(п) = ХА се(ј)= 2 Ave(n-i) 
2-1 1=0 
(4.8a) 
and 
nae O y(n) + e'Cn) 
n i 
- C' E A'e(j) * e'(n) 
j=l 
n=l , 
= C! 2 дае (12) + e' (n) CH 8p) 
1=0 


Equation (4.8a) is identical to Equations (4.6a-c) but 
(4.8b) has a subtle difference. This difference is that, for 
the T forms, the output error depends on all previously 
generated errors in the state equations including the most 
recent (i.e. at time nT). The error Av(n) still only depends 


On the most recent error e'(n) of the output equation. 
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miemorrrerence in the expressions for Ау(п) for the 5 
and T forms is a result of the choice of state and output 
equations and will present no difficulty in the stochastic 


development which follows, 


С. STOCHASTIC FORMULATION WITH CORRELATION 

A general expression for the variance of the output error 
when the quantization errors are assumed to be discrete ran- 
dom processes can be developed by making use of the error 
propagation equations of the previous section. This expres- 
mien сап then be used for analyzing the specific cases of 
the second order filter and particularly the canonic forms 
pr Chapter III. 

1. General Formulation 

Consider first the T forms of the state equations and 

thus Equations (4.8) as expressions for the accumulated 
Abate and output errors. Assuming for the present that all 
Random variables and processes are zero mean the variance of 


mme output error can be expressed as 
2 E L y 
o,y(n) = ЕШШ (т) Л ГАУ (ПОК (4.9) 
where E[:] is the expected value. 


P ppl ying the first equality of Equation (46860), 


the output error variance becomes 
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слуп) = ЕГ(С'у(п) + є'(п)) ((є'(п))° + уп) (6')%)] 
= C'A (n) (C')Ë + С' ELy(n)(e"(n))*} 
+ Ее! (п)у (п) 1( с") + Ее! (п) (е! (п)) 81 
< СА (а) (с!) + 2с" ЕГу(а)е! (п)1 + o, (п) 
(4.10) 
where 
А (а) + ЕГу(а)у (а)! (4.108) 
с (п) = E[Ce! (n))*] for the single output case 
(4 105 
апа 
ЕГе! (п)у (а) (с!) + с"ЕГу(а)е! (п) 1 (4,100) 


Equation (4.10c) arises when a scalar output is 
assumed; thus the cross terms in (4.10) are combined in the 
final expression of (4.10). 

A (a) is the covariance matrix of the accumulated 
state errors at time nT. The dependence of К on the 
Meecistics of the quantization errors is found by substituting 


ENS) into (4.1038). Then 
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п 
5 


й 
Е А-у ео) 27305 
y = 4-1 
= Gey ee t.nej 
= E ; А "Ele(i)e’(j)](A°) CUSTI) 
1=1 j=l 


Defining a quantization error correlation matrix as 


R (n,m) $ Ele(n)e"(m)] (4.12) 


Prepiles (4.11) to be written as 


n n n-i tones 
MA AR AAA Е 
y 1=1 ігі е 


Tt is now possible to write the expression 


n n E 
— е . S zb 
151 11 


) E 
т 


n-l n-1 : MS 
РТУ AA (ТОКЫ, J (4.14) 
i=0 1-0 


Subtracting (4.14) from (4.13) yields 


Ca ne Can 
A, (n) e АЛ (п)А = R, (n,n) - А R,(1,1)(A ) (4.15) 
n-1 и: n-1 s 
г В (п,ј) (А0) APRA ma) 
5 е А е 
Jal Je 
n-1 n-1 i 
A даа in (ана (до 
1-1 = 1-21 
n-1 n-1 tn 
H z АТН (4,3) = ReO (47) 
ісі j=l 


Tol 





Equation (4.15) can be simplified if it is assumed 
that all the processes involving quantization errors are wide 
sense stationary and mutually wide sense stationary. Under 
Bais assumption, the correlation matrix Re (n,m) is a function 
only of the difference of its arguments, (m-n)(taken for 
notation purposes in that order). Then the correlation 


matrix can be written 
R (n,m) = R, (m-n) (4.16) 


Under the assumption of wide sense stationarity, 
the first and second terms of (4.15) contain R,(0). Further- 
more, the double summation term becomes the 0 matrix, so that 


(4,15) reduces to 


A (n) - АА (од: в_ (0) - АПБ (0) (Аб)? 


n-1l4 а 
е, R. (J-n) (A)? 7 


jai | 


APR (3) (a0) 27 


5 
iba 
HE 


1 


n-l ^ 
i APT R (-1)(A 


1 
CH 15 


TER (n-i) _ E 
al b 


Now consider the entries of the correlation matrix: 


eh and components of e(n) 


ana e(m) respectively then the i entry of the correlation 


Lf e, (n) and e, (m) are the i 


Matrix is 
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[R (n,m)];; = [Ele(me (m7), (4.18) 


Ј 


= e, (2e, (m) - e, (m)e, (n) 


t Е 
ГЕГе(т)е (п) 17,, = (8, (т,п)2,; 


This implies that the transpose of the correlation matrix 
is the correlation matrix with its arguments reversed. For 
the case of mutually wide sense stationary processes this 
means 


^ 


ү 
Ке(к) = Ке СК) (4.19) 


Making use of (4.19) in (4.17) it can be written 


bnat 


t © ja tren 
А (п) - АЛ (п)А R¿(0) - А R,(O)CA ) 0220) 


п- 1 А А 
+ X {ATTUR, (n-i) + VS OA 
i=] i 


n-1 oA ^ 
x = : С - 
- © {APTĪR (na) (ah)? + [APT R Ea a 
i=l 
Now as n becomes large the last summation of (4.20) 
becomes negligible since the filter is assumed to be stable 
which implies that lim A” = 0. The second term in (4.20) 


n> 


also becomes negligible so that in the limit 
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> ща 





A, = An A? = R (0) + E (АТ G) n CatR, (1) 1°) (КОО, 
pr 


(ATR, (1) + R (1) (A (1.215) 


it 


и са 8 


- АЛ А 0) + 
Ay J Re! ) 


1 
The left hand side of Equations (4.21) can be solved 

for AY so that these equations essentially give an expression 

Бог Аз the covariance matrix of state errors, in terms of 

the correlation matrices of the quantization error sources 

for large n under the assumptions of zero-mean, wide sense 


stationary and mutually wide sense stationary error processes. 


t 
- ANA’ = Е 4,22 
A, А, ( ) 


(where Е 15 a symmetric matrix which represents the right 
hand side of (4.21)) can then be substituted into (4.10) as 
Mert Of the solution for the output error variance. 

The factor E[y(n)e'(n)] which also appears in (4.10) 
must be found before the output error variance can be written 
merely in terms of the statistics of the quantization error 
Sources. From (4.8a) it follows that 

Inc 


Ely(n)e'(n)] = 2 A'Ele(n-i)e'(n)] (1.23) 
1=0 





Now define a cross correlation matrix (in this case 


a vector) between e(n) and e'(m) as: 


R, Спут) 8 ЕГе(п)е! (т) 1 (4.24) 

Under the assumption of wide sense stationarity, as 
More, the cross correlation vector depends only on the 
difference in its arguments, 1.е. 


^ 


R, i (n, m) = Ro ¿1 Mon ) (4,25) 


so that (4.23) becomes 


Ely(n)e'(n)]= Е АВ, , (4) (4.26) 
u ЩЕ қ 
End it follows that 
ШЕУ (0) е (01 = x MR, (09) (4.27) 
n> 1+0 ~ 


Substituting (4.27) and the solution to (4.22) into 


ШТО) yields the output error variance for the T forms, 


jt T РН. 


E EN 
E 


R O 
120 55 


2 1 = 1 1 
ШІ 6.6!) < С AS 


As a function of the statistics of e and є'. 
For the S forms, the difference in the output error 


propagation equation (4.7c) causes a slight cifference in the 


expression for the output error varlance. 


Що 





Substituting (4.7c) into (4.9) yields 


2 
Sayin) 


Е[(Су(п-1) + є(п))(є (п) + у (п-1)с°)] 


2 
E 


CA, (n-1)C* * 2C E[y(n-1)6(n)] + о (4.29) 


The first term of (4.29) contains А (а-1) instead of А (а) 
Which appeared in (5.10). But in the limit 

B (n-1) - Л (а) = ло The second term of (4.29), however, 
results ina slightly different expression in the summation 
form. Using (4.7b) 

2 


г дЇЕГе(п-1-1)є(п)] 
1=0 


n 


Е[у(п-1)=(п)] 


n-2 1^ | l 
= 2 А Ro (1+1) (1.30) 
1=0 


which in the limit becomes 


e ^ 


aL | 
КЕНЕ (4.31) 


И 68 


ПИ = Јел) | = 


no i-0 


Ehen the expression for the output error variance for the 
5 forms is 
2 


2 t T . 
--. + 
с, (е,Е) = CA (ЛС е оС А Ке (+ 1) + 0. 


(4732) 


Equations (4.28) and (4.32) together with (4.21) are 


the principal results of this section. They express the output 
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Ж ИЕ Variance of the S and T forms in terms of the statistics 
mt the internally generated quantization errors including the 
effects of correlation. The assumptions made were that all 
processes had zero mean and were wide sense stationary and 
mutually wide sense stationary. 

When the assumption of zero mean processes is relaxed 
the expressions do not change appreciably. All that is 
necessary is to redefine the covariance and correlation 


matrices: 


А (а) А Е[{у(п) - my) Hy(n) – m(y)}*] 
= ELy(n)y"(n)] - m(y(n))m"(y(n)) 
= A (n) - m(y(n))m” (y(n)) (4.33) 
#_(п,п) 8 ЕЦе(п) - п(е(п))}{е(т) - mtetm))}] 
- ЕГе(п)е (п) 1 - m(e(n))m"(e(m)) 
= R (n,m) - m(e(n))m"(e(m)) (4.34) 
R (n,m) ЕЕ = m(e(n))}{e(m) – и(е(т)))) 


| 


E[e(n)e(m)] - m(e(n))u(e(m)) 


H 


R,¿(n,m) - m(e(n))u(e(m)) (4.35) 


T 





where m and p are the means of their arguments. 
When it is assumed that the processes are vide sense 


stationary it is possible to write 


A (e) = А (е) ~ т(у)т (у) (4.36) 
R (n,m) < Қ.(пеп) « Қ.(п-п) - п(е)п е) (4.37) 
pum R..G-n) = R_ (men) - m(eJu(e) (38) 


where the dependence on a time index has been removed because 
wide sense stationarity implies dependence.only on time 
 ШТТІегепсев. 

Then the expressions for the output error varience 


for the S forms is 


2 Қ” т ам 2 я 
О ду(е,Е) < СЛ (ес + an Ro zen 64 550) 


and for the T forms 


2 F е o 2 
g av езе") = СА, (е) (С!) Ee | “| Reet Ci) + Go, (4.40) 


al 


where Ay Ce) "пе solution to 


бои но > 
; (А R (4) БИК CON ) 


- Ah Пар (о) + 
Ay(e) - AA, Ce) (0) +. 


1. 
(4.41) 
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апа 


2 


o_ = El(e - и): = ЕГе 1 Um 


(4.42) 

These equations ((4.39) through но леге 
general expression for the output error variance assuming 
Endy Wide sense stationarity. The effects of structure are 
imbedded in the nature of R (1), CY and с and in the 
Entries of the matrices of the state and output equations. 
The next section examines these effects. 

a. Structural Effects (General) 

In order to take topological or structural effects 

into account, it is necessary to consider the terms of the 


Borreiation matrices cf error sources as given by (4.16) and 


Ц л 
в. 25). When zero correlation is assumed, В. (К) 15 Бакеп 


to be zero for all k # 0, and Е. (0) is taken to be a diagonal 
matrix where the ва entry on the diagonal is the sum of the 
second moments (variance) of each individual error source 
comprising e, (n). Assuming uniform distributions, these 
terms become? (0°) times the number of unique non-zero non- 


th row of the matrices which 


unity multipliers in the k 
contribute to the error e, (n). 
When correlation is assumed (both spatial and 


time) the picture becomes much more complicated. In general, 


“This development returns to the assumption of zero mean 
processes for ease of notation since there is no loss of 
EEnerality in Going so. 


2g* = һ2/12 for rounding, h^/3 for trúncation. 





the (1j)*h term of the correlation matrix is glven by 


СОЯ = [temer] = ELe, (n)e , (m) 
р; а, | 
= E A масы ыы т 
р. 9 
= dí qe Ley inde y (m) ] (4,543) 


мһеге p, is the number of multiplier errors concr rIDUtCINE TECO 


e а. 15 the number of multiplier errors contr outing to 


1? 2 
E and k,£ are additional indices to identify the contribu- 


Erons to е.» е; respectively: Пар S 


р. а. 
i y 
е; (п) = i е. (1) апа e, (n) - м (4,44) 


In all, then, there are six indices to be considered to 
metermine the entries of the covariance matrix; two (i,j) 
associated with states; two (k,%2) associated with multipliers; 
and two (n,m) associated with time. When wide sense station- 
arity is assumed, then one index representing time difference 
peplaces the two indices of time. 

Considering a typical contribution to Е (пт), 
such 25 Ele, (n)e, y (tm) ], the question to be asked is when is 
this correlation statistic non-zero? One obvious answer is, 
Ehen the multiplier coefficients are the same number and the 


Signal quantities are identical (e.g., x» (n) - x, (n-1)). Hera 
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о case there is a unity correlation coefficient and 
E 2 2 
ELe (m)e,y(m) ] = Ele ¡y (Mm) JELe ¿y (m) ] : (4,45) 


Another possibility to be investigated is when 
mee Signal quantities are identical, but the multiplier 
coefficients are different. If the signal quantities have 
 (55шап distributions then the correlation coefficient is 
a function of the ratio of the two coefficients as shown in 
Appendix В . AS expected in this case if the multiplier 
НИ cientus are the same, the error statistics have unity 
correlation. Note that the correlation coefficient depends 
пе Sign of the ratio of the multipliers. 

When the multiplier coefficients are the same 
But the signal quantities are not identical, then it is 
Suggested that the correlation of errors will depend mainly 
on the correlation between the signal quantities. 

Finally, the most complicated case is when the 
moetrficients are different and the signal quantities are 
Mmereidentical. It is expected in this case that the 
Amr elation is small. 

The foregoing discussion applies equally well 
to the terms comprising ет). Experimental determina- 
tion of correlation effects are considered in Appendix D 
for second order canonic filters. In the next section the 
general error variance equations are simplified for the 


second order filter. 
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2. Simplification tc the General Second Order Section 
a. For Arbitrary State Equations 
For the second order filter the solution to (4.22) 


is a symmetric matrix 


À Х 
Уі1 Y12 
ци Е (4.46) 
y А X 
where 
Е 2 2 
br. г. са 811855) (1-855 ) = 8, 5,854 (1*a5, ) 
аа па а 5 (4 46a) 
22551292052 и 22 
+ £__[(lta..a аа я = 
22 11550 та аа 
\ Е) — (ЕРЕ | а “уа а + а „а “a ] 
y y id ee 2] 122 2 
то 2 
2 Е 2 
+ £, „С1-а 1) (1-а,5) -- 8.5854] 
* f..[(1-a SR a + а а Ca TIAA още еди 
22 77% 12 00 Mss oc oT 
ЖЕ еве аа зарына. айар “1 
у (a 11°22 Толе Би 
а аа ба а (а) (4.460) 
Ma 112 2] A 11 


2 2 
+ f55L (1-844855) (1-844) - а. „а (1ta,,) 11/4 
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with 
А = (1-а IM ы а - а.,а,.) 
ШІП! 22 o le al 
2 2 
2(1-841)84 5851855 
E En. > а 
2(1-855)841845851 - (1-84,85,-845851)8,5857 
2 2 
Е Ча уа а ,3а.., (4,464) 
and 
“Е 
Е = . . with fio = fo] (4. 46e) 
2) 22 


For general second order state equations this 

solution could then be substituted into (4.28) or (4.32). 
Оши ог Canonic Arrays 

When the state equations are specifically those 
(Пе canonic arrays in Chapter III, it is possible to 
Eubstitute the coefficients found there into Equations (4.46) 
Eo find E in terms of those coefficients and the entries of 
К благах в. If the factor Ely(n-1)e(n)] in (4,29) is 
Biven by 

E . 


a= ЕЕ е ра E 209 (4,47) 
Es г 10 
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then the expression Гог the output error variance for the S 


forms of the state equations become: 


= 


SM: c = (file + c + b) + ce(-2a)] 


+2 f Pie + SED t се Па Ел 


12 


+ го бе Бе (а+) - ?ce(db^) + ес рблев) – а карате 


2 
+ 2(св. + ев.) + А (4,48a) 
е 2 
е = + ыы 
SM”: B, (£4. 0 b) да > + До СЕЗ (4.48b) 
SN: Au dM * e^) (14b) + co(-2a)] 


+ des т e“)(1+b+ab) + се(1-а°-Ь°-2а) ] 


+ г „Се (21+ъчаь) - (145) (1-ь)) 


* 2ce((1-b)(14b*ab)-a(1*a*b)) 
Тт о. 200 = нео (про | Пи 


* 2(c(g, * Bj) + е(62)} + 0% (4.480) 


| 


оо p 
SN': c {22.1 (1+а+5) 2Р. >(1+а+Ъ) + Р52(1 b)378 


(4.460) 
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where 
A = (1-b) [(1+b)° -. E (Ups 


Mnen the T forms of the state equations are used, 
ШІ (с factor Ely(n)e'(n)] in (4.10) be given by 
51° 


С! = БЕ) е mo 
Во. i 


и с 8 


SM 

A Ree Фекер» 
0 
For the T forms, the expressions for the output error variance 
are identical to Equations (4.48) for the coresponding canonic 


arrays (e.g., TM corresponds to SM, ТМ E to SM | 


ессе 
ос that c, = Е.» Во апа of are DUC Cd DY нем Ec M 517, 
B5' апа is Pespectlvely. 

Note that the transpose forms do not apparently 
епа оп с апа е. The factors E however do depend on 
these coefficients since in the transpose case they appear 
in the state equations rather than the output equation. 

The quantization errors they cause will appear in e(n) thus 
in ‹ and thus in TTE 

The output error variance for the transpose cases 
does not however depend on any G or of since e(n) (or & (n)) 
ШИ (сто for the transpose case, there being no non-zero, 
non-unity multipliers in the output equation. 

The expressions for the output error variance 


given in Equations (5.58) are identified without subseripts 


but rather by an identification which implies the applicability 
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of that expression.to all members of the set of realizations 


Eh that basic identification (e.g. SN is a member of 


10 
the set of SN realizations). The reason these expressions 
apply to the entire set is that the development to this point 
has depended solely on the basic state equations. The dif- 
ferences of structure implied by the difference between the 
eight basic forms of the state and output equations exhibit 
themselves in the four basic expressions of (4.48) for the 
S forms and the same four expressions for the T forms. The 
differences of structure imposed by the various manipulations 
of the flow diagrams are entirely embodied in the factors 
E and бі. What remains to be examined is how that structure 
16105 itself in those factors. Structural effects were 
Giscussed in general in Section C2. Now they will be examined 
in the context of the canonic arrays starting with the simpli- 
Bang assumption of uncorrelated error sources. This condicion 
и111 then be relaxed to achieve some better insight into the 
EDfects of correlation. 
пиши ucorrelated Case 
Wnen all the error sources are assumed to be 
uncorrelated, the right hand side of Equation (4.22) becomes 
Eles(n)] 0 
aq = (4.50) 


0 Е[е5(п)] 
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2 
where Е[е; (п) ] is the sum of the variances of the individual 


multiplier errors contributing to e,(n). 
For example, in the realization of the array, 


B. the first state equation is 


x4 (n) = с(х. (п-1) * u(n-1)) - (atc) x, (n-1) + X5(n-1) 


(4.51) 
therefore 
e4 (n) = etx#(n-1) Ate = [cixf(n-1) + u(n-1)}], 
+ {-(а+с)х* (и-1)} -. [{-(а+е)х1 (1-1) 1. 
% х5(п-1) - [x3(n-1) J, 
= e (n) + е сас) (2) + 0 (4.52) 
and 
ЕГе (а) - ЕЦе (а) + е сане) (9011 
2 2 
= Ele, (n)] + Ele a4.) m] (4.53) 


since Ele (пје y(n) J is assumed to be zero because the 


(ate 


error sources are assumed to be uncorrelated. 


167 





If one assumes that all error sources are 
identically distributed with variance сё, then for the 


example SM the variance of e (n) 18 


20 


Elef(n)] = 20° (4.54) 
In general, for the uncorrelated case 


B (4.55) 


where Ще (1 = 1,2) is the number of unique non-zero or 


th 


non-unity multipliers in the i row of the canonic array. 


For the example, the array is 


КО г; 1 с 
E 
SM e E Ер 0 е (4.560) 
ШІ 0 а 


It appears at first that there are three non-zero, 
non-unity multipliers in the first row. But by referring 
to the flow diagram of this realization (Figure 3-20) or 
Equation (4.51), it is seen that there is only one multiplier 
ОП с as the coefficient and one with atc. Thus there are 
two quantizations required, one after each of these two 


multiplications. 
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In the array 


с-(а+с) -b Ш 
Mo = 1 0 0 (4.57) 
с е а 


Џлеге аге three unique non-zero, non-unity coefficients in 


the first row and none in the second. Then for SM 


20 
30° 0 
рле (4.58) 
0 0 
whereas for SM t 
20 
20° 0 
P s 
0 өре” 


meeretore in counting coefficients for forming F, only 
unique non-zero, non-unity coefficients in each.row are 
included, because the same reasoning is true in all the 
ша! 1га 1015. 

То find the entries of the vector G or G', it 
is necessary to examine the relationships between the state 
equations and the output equation. 


From Equation (4.47) and the definition of Ruck), 


G = Z ME[e(n-1-1)e(n)] (4.60) 
150 
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and from (4.49) 


AlEle(n-i)e'(n)] (4.61) 


о) 

| 
и са 8 
e 


At first glance, Equations (4.60) and (4.61) 
would appear to be zero when error sources are assumed to 
be uncorrelated. But it is to be remembered that e(n) and 
e(n) are not error Sources. Rather they are composite errors 
made up of the errors from possibly several sources. In 
fact, in some of the realizations of Chapter III a particular 
Arce of error may find its way into a component of e(n) 
and also into e(n). 

Consider for example the flow diagram for tne 
pcalization TM. (see Figures 3-20 ). The first state eavation 


д 


mor this realization is 
х (п) = -ax, (n-1) - bx, (n-1) + u(n=1) (476242) 


while the output equation is 


v(n) c' x (n) = ах. (п) + (ate! )x4(n) вао 


n 


c' x4 (n) = ax, (n-1) + (ае! )хо(п) + d'u(n-1) 


(4.625) 


Тһе еггог e, (n) in this case represents the sum 


of the two errors from Guantizing ex, (n-1) and bx, (n-1) while 
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е'(п) 15 made up of the errors from quantizing с!х. (п), 


(ate')x,(n) and ax, (n-1), so that 
e, (n) = e,(n) + e, (n) (4.63a) 


and 


е! (п) е 1 (1) + е (а+е!) (п) + e, (n) (4.63b) 


Then 


Ele, (n)e'(n)] = Ele¿(n)e,,(1)] + Ele, (m)e ¿401 ) (097 


+ E[e*(n)] + Ble, (nje,, (ns J 


+ Ble, (ne 00] + Ele, (n)e, (n) ] 


(4 .64a) 


(ate 


Under the assumption of uncorrelated error sources all these 


terms are zero except Ele“ (n) J so that 
Ele] (n)e'(n)] = Elef(n)] # 0 (4.646) 


ES implies that in Equation (5.61) the first term of the 
EUumation is not zero. The other terms are zero, however, 
in the uncorrelated case. 

Пише хтео ог а1] realizations it turns 
апа ТМ 


out that only ТМ ras pattern. А11 transpose 


|. 2 
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configurations require that e(n)(or e'(n)) be zero since 
Mere are no non-zero, non-unity multiplications in the output 
equation. For the SM and SN forms, although in most cases 
there are coefficients which contribute to both e(n) and 
E(n), the term E[e(n)e£(n)] does not appear in Equation (4.60). 
TMo апа ТЫ have no coefficients which contribute to both 
men) and Е' (п). 

Again assuming identically distributed error 


^ 


sources with variance о. then for the realizations TM, and 


ik 
TM, 
f ^2 
51 0 
JE | = (4.65) 
52 0 


Although there are only two cases when G! is 
non-zero for uncorrelated error sources, the exercise in 
finding them will be helpful in understanding the effect of 
meructure when correlation is considered. 

Finally, the last term of (4.28) or (4.32) 
involves the variance of the quantization errors contributed 

2 2 y. 


By the output equation alone, 1.е., Oc (or Tot 


When 
uncorrelated error sources are assumed, this term is the 
В Of the variances of the individual errors. Assuming 


Meaentically distributed error sources, 


с“ ES - (1.65а) 


2 





where v is the number of non-zero, non-unity elements in 
the output row of the array. 
Using the previous results, Equations (4.48) for 


the S forms become: 


SM,SN: ee = {(u,/d) (еже?) (1+) - 2ace] + vo? 
(4.662) 
сме: Ja = (и, Е M4) (1+b)0°/A (1.665) 
ЗМ: с, = {2р. (1+а+ь) + и (1+6) } 02 ил (4.66е) 
5 Ау 1 2 қ 


where A is defined in (4.48e). 


For the T forms these equations become: 


ТМ ТИ: О, = ИШ ОС i (е')2) (1+) -2асте! ]+0)02 
(4.6643) 
TM. ,TM,_: of = И Иво ги Сен None) со a ее 
| 2" Ау 1; 
(4. ббе) 
ШЫ a u a a bo AA (1.667) 
р Лу 1 2 р 
ТМ": c? = {әр (1+а+ь) + и. (1+6) }02/А (4.666) 
j Av al 2 7 


ые гаст that H = O for non-transpose configura- 
tions is already included in Equations (4.66). 

Maple 4=1 lists the results of the analysis 
above for each canonic array Гор theluncorro ated m identically 


E Stributed case. 
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АВКАУ a aie: ӨР» 

А 22 

SMo SN oo Due (2p+2)0 1 
sa 

SMo >5M, y SNo1 57 2009 3 (2p+3)0 2 
да 

SMo5,,SM5g,SNo5,SN,9 3 0 2 (3p+2)0 3 
^5 

SM, 5 ,8M5, ,SN4 5, SN 54 3 0 3 (3p+3)0 4 

20 _ 

5М,1>5М44 2 0 Ц (2p*l)c 5 
= по 

SM, ,SN,,, . Ц 0 2 (Цр+2)а 6 
uU Na 

SM, , TM," 2D е ц (1+) с /А 7 

sn... my? 2 2 D [6(14b)4lla ]o^/A 8 

4j? ag 

TM, TN. 2 0 2 (2p'+2)0° 9. 

TM, 2 Q0 3 (2p'+2c!+3)0” 10 

TM, 3 0 2 (3p '+20!'+2)0° 11 

E 2 2 
р = [(ес +е )(1 + b) - 2ace]/A 
t = rye rye tat T/A 
mes (ет) + (e')~)(1 + b) 2 2ac'e* j7 
Men = pyr b) = acd 
TABLE 4-1 2 FOR UNCORRELATED, IDENTICALLY 


FORMULAE FOR o 5 


^2 
DISTRIBUTED, ERROR SOURCES (MEAN=0;¡VARIANCE=0 ) 
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Appendix С verifies that none of the formulae 
for the uncorrelated case violate the requirement that they 
predict a non-negative number for cee Briefly, the results 


of that appendix are that for a stable filter 


А > 0 (4.67a) 
peor p') > 0 (4.675) 
ПО > ~ i (4.67c) 
(1+b)/A > 1 (4.674) 


Examination of the relative magnitudes of the 


Expressions for с 2 shows that the smallest value for formu- 


Av 
ESe 1 through 6 is (2p + 2ус2 and from (4.67b) the smallest 
this can be is 20°. Formula 7 is at least 102 busco 
and formula 8 is at least 202 since a is at least -(1+b). 
Formulae 9 through 11 predict a minimum variance of 20°, 
2.502 and 1.502 respectively. 

It has been stated empirically | J that wor 
the canonic realizations derived by Parker and Hess, noise 
generation is less in the S, (SN) forms than in S (SM), but 
that the transpose forms had still less, S/(SN") being 
better than S^ (SM), i.e. the least error occurs in Sn”, 


Mable 4-1 does not predict this, пог does it contradict it. 


But formula 1 predicts that SMoo and оо yield tne Томе ще 


ds 





noise for the untransposed S forms, making no distinction 
between SM and SN. Formula 8 says any of the SN" configurations 


are probably as good as SMoo and SN... Furthermore, that 


00 

same formula indicates that the TN" forms are comparable 
and formula 9 adds TMo and IN. to the group with a minimum 
noise of 202. 

The most intriguing suggestion of Table 4-1 is 
that IM, may yield the lowest noise, having a minimum variance 
БГ Ш 502. 

But all this is highly conjectural especially 
since only minimum values are being compared. There are 
ШЕШ sets of coefficients which will result in these minima 
being achieved. For example, (1+b)/A is only minimum at 
КОШЫ кїп in the (a,b)-plane. The factors p and p' cam be 
minimum for any value of a and b as long as c and e (ор с! 
апа е!) are both zero which is a degenerate case. 

Furthermore, this analysis was based on the 
Eumption of identically distributed, uncorrelated error 
sources. More useful and interesting results can be achieved 
by relaxing these conditions. Furthermore it is closer to 
reality as will be seen in the next section. 

cor related Case 
When the quantization error sources are assumed 


to be correlated, the quantities Г.) in Equations (4.48) 


are given by the matrix equation 
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Ес В. (0) + : 


(a*R (4) + LAT R (29) (4.68) 
i 


1 

The entries of RJ were discussed in general 
in Section C-1-a of this chapter. It was shown there that 
these entries consisted of auto- and crosscorrelation terms 


of the error sources. 


For example, consider the realization тм Lor 


which the array is 


-а 1 с! 
-р 0 е! 
16 0 da 


and the structure is that or Figure 3-42. The error vectors 


mor this realization are given by. 


e, (n) e,(n) + ес! (п) 
= (4.69a) 

егіп) ey (nj + e,,(n) 
е(п) = 0 (4,695) 


where е „(п) is the error introduced by the multiplication 
of x4 (n-1) by à and the subsequent quantization of that pro- 
Euct at time nT. The errors ey (n), ес! (п) апа есі (п) аге 
similarly defined. Figure 4-la shows a schematic represen- 


mation of the error sources for И, Figure !-1b shows the 


ү 





composite errors e, (n) and е> (п) as composite error sources 
where e (n) end e 5(n) are given by (4.69a). 


The correlation matrix R, (1) is given by 


Ele, (n-1)e, (n)] Ele, (n-1)e,(n)] 


R n» = ELe (n-1)e" (n)] = 
а E[e,(n-1)e, (1) Ше 1)e, ma 


(MERO ) 
westi tuting (4.69) into (1,70) УД61045 
a 8114) В. (1) 
R, (i) = (1.72) 
Ro (1) 8-01) 
with 
Кес) : Ble, (n-i)e,(n)] | (4. 71a) 


(1) + Е (1) qub) 


END - ROG) + R ,(1) +R m 


HT aa ас! 


о-в (i) + Е (1) + В (1) + R (15) (Да те 


ща ар ае! ср сте! 
Roy (1) = E, + А ос! (3) + Rota lh) + Ке о: 1) (4.714) 
ВЕ = Rp (1) + БЕСІ) а Rot, it) + Rote (i) (4,71e) 


where 
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Rag 3) 8 Ele¿(n-1)e¿(n)] (4.71f) 


under the assumption of wide sense stationarity. 


Appendix D describes an experiment which was 
conducted to check the validity of Equations (4.48) for the 
realization ШТ. The "actual" output error variance for a 
white noise input was calculated by summing the square of 
Ene difference between the outputs of a double-preeision 
Е.а and one which rounded products во опе deem 
ШӘП. (This, of course, required the assumption that the 
process was ergodic (or stationary in a strict) sense in order 
for the time average to be equal to the ensemble average.) 
Baring the process of rounding products, the error from each 
Quantization was saved, correlations were calculated and 
substituted into R. (1) and F was computed. The "predicted" 
output error variance was found using (4.48b) and the a 
posteriori information about the correlation of the error 
MT ces. Excellent correspondence between the actual and 
predicted values was observed for three sets of transfer 
function coefficients. 

Of course F was not computed as an infinite sun, 
but a sufficient number of terms were used to observe the 
convergence. Of particular interest is the observation that 
me higher damping coefficient caused more rapid convergence. 
This is to be expected since it is the nature of a stable 


mercer that as the power of the matrix A becomes larger the 


179 





entries converge to zero and that the rate of this conver- 
gence depends on the damping coefficient. But each term of 
F also includes a correlation matrix. The experiment shows 
that convergence to within 0.1% is faster when the filter 
is driven by an uncorrelated input than when driven by a 
highly correlated one such as a step (Heaviside) function. 
This too is to be expected since a highly correlated input 
implies higher correlation in the states and thus higher 
ErPrelation in the error sources. Tres tends to make) tne 
Entries of the correlation matrices larger, hence a slower 
convergence of F is expected. An heuristic examination 

Mi this idea follows. 

Consider the nature of the auantization error 
source. Figure 4~2a shows the function find), ПОН ЕС 
where f(n) is the number to be rounded and h is the quanti- 
Eation step size. The dotted line represents the identity 
Bunetion f(n) » f(n). Figure 4.2b shows the error 


 Т(п) | = f(n) - [end]. Note that the error can be written 


as 


|| 


е(п) = е(#(п)) = f(n) - hk(n) (4.72) 


with 


k(n) » [(£(n) * 2/hJ, (4.73) 


where [+], means the largest integer less than (+). Equation 


180 





(4,73) implies that 

(f) - 2/h. € k(n) « (fin) + Пу (4.78) 
which implies that 

bk(n) - 2 « f(n) < hk(n) + È (4,75) 


Now consider the autocorrelation of the error 


E[e(n-i)e(n)] * E[(f(n-1)-hnk(n-i))(f(n)-hk(n))] 


E[f(n-i)f(n)] - hnE[k(n-i)f(n) * f(n-i)k(n)] 


PESE a 


E[f(n-i)f(n)] - hE[f(n-i)k(n)] 


- h(E[k(n-i)f(n)] - hE[k(n-i)k(n)]) (4.76) 


If the discrete probability density function 
пра. г.) of the random variable f(n) is given by pp(f(n)) = 
ШО = f) and the p.d.f. of k(n) is py (k(n)) = РШ x 


then by (4.74) and (4.75) 


py(k(n)) = Prf(f(n) - 8)/n < k(n) < (£(m) + 3)/n] 


| 


һ h е 
P [hk(n) - 5 < fin) < hk(n) + 5) (4477) 
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If the cumulative distribution function (c.d.f.) 


ШІ f(n) is given by Pp (GF Q)) then 


Рк(К(п)) = Р-(а,) - Рь(а-) (48) 
where 
Y h 
x ELE = (4.78а) 
B h 
as = hk(n) + > (178b) 


Шоле says that k(n) is distributed in such a way that its 
Statistics are very similar to those of f(n) so that Equation 
(4,76) looks like it has four terms that are very similar 
ша form to E[f(n-i)f(n)], the autocorrelation of fín). That 
is, it appears that the form of the autocorrelation of the 
ОГ should be similar to that of the autocorrelation of 
шие" пипрег Со Бе quantized. 

In Appendix F it is shown that the autocorrela- 
Nn of the state variables of a linear filter (i.e., an 
infinite precision linear filter) can be found, given the 
Mitocorrelation of the input and the transfer function between 
the input and the te er consideration. In the previous 
Paragraph it is suggested that the autocorrelation of the 
errors is related very closely to the autocorrelation of the 
ENDS] guantity to be quantized. In all cases this signal 


Entity is the product of a coefficient and either a state 
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me une input. The process of multiplying by a coefficient 
Mees not effect the autocorrelation of the signal quantity 
before or after the multiplier except to scale it by the 
 Паге оі the coefficient. Therefore, it is suggested that 
Æe autocorrelation of each error is similar in form to the 
autocorrelation of either a state variable or the input. 

A similar argument suggests that the cross correlation of 
Гог sources is related to the cross correlation of two 
signal quantities (or to the autocorrelation of a single 
signal quantity if the two coefficients multiply the same 
Msigenal quantity). 

For simplicity, two extreme cases that bound the 
forms of correlation functions will be considered: (1) when 
the input is uncorrelated and (2) when the input is highly 
correlated (i.e., its autocorrelation is a constant). 

(1) Reduction for Uncorrelated Input. Appendix 
F shows that when the input is uncorrelated, the auto- 
correlation of the states will be a function bounded in 


magnitude by ge" YInl7 


where K is the variance of the state 
and y is proportional to the damping сове ел. ов обе 
transfer function between the input and the state under 
consideration. But by Equations (3.9) and (3.11) it can be 
seen that the characteristic function for the state equations 
is the same as that of the overall transfer function so that 
Y is the same for both states and output. 


The cross correlation between the uncorre- 


lated input and each state is proportional to the unit pulse 
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response of the respective state and the cross correlation 
between the two states is proportional to the correlation of 
the two unit pulse responses. All of these functions are 
then bounded in magnitude by an envelope proportional to 
erialt 

Assuming the correlations of the error 
sources behave in a manner proportional to the signals before 


the multipliers, these correlations may be written as 
ag 
Ав (221 < Кове (4.79) 


where Rag 3) is defined as in (4.71f) and Kag is some positive 
proportionality constant involving the variance or covariance 
ПИ errors and other factors. 

ТБ may be possible to substitute (4.79) into 
a form of Equations (4.41) which involves the magnitudes of 
the factors f,, and g,. “йеп сотая ЈАН апа le, | one 


could examine 


ВВК (0)> + Z {<a*><R_(i)> + катък (1) 16) (4.80a) 
у ЕШ 
Ex E «AT »«R (1+1) > (4.80b) 
— ев 
із0 
En!» < I <АТ><В ЕР (4.800) 
E ес 


where < > indicates taking the absolute value of each entry 


of the matrix contained therein. This process may however 


18! 





give away tco much to obtain a good estimate of the output 
error variance. This avenue is not pursued here but is 
suggested for possible future study. 

(2) Reduction for Highly Correlated Input. 
Appendix F shows that when the input is highly correlated, 
the autocorrelation of the states will also be constant. 
Assuming the error sources also have constant autocorrelation 
and cross correlation, the matrices Ri), Наси (3) and 


Ке (1) in (4.21), (4.28) and (4.32) will be constant for all 


mee Рог this case, let 


Q Q 
Қ Sis ig 
R, (i) = 92е = (12812) 

Q Q 

|2251 SON 

R Е í 1.8 
В (1) Qog! Ей Че! IE ( e 159, 
B Ges S 5 (4.81с) 
Ке (4) B des ü "Че e jd 


Then the matrix F whose entries appear in 


(4.48) can be written as 


NE A e 
F=Q + 2A Q en Q, J 
6 1451 4-1 
2: t 
= 0, + АГ1-АЈ ТО, + (GU-A178,) (4.82) 
By the matrix identity 2 А = ТУА. Тһеп 
n=0 
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( аа” 
Qe e, (1*8411785278,,855*8,585,) * 28,4Q, , | 
Ж JU zen) 
fuc (4.834) 
l tat b 
(0 ео уа © Фа 0 +(ал-а. а,„а „ад. )9 ) | 
Е _ е е» 22 12 ее: 21 ее, РЕД "а ага Оно ее, | 
ES ЕЕ | 
(1.835) | 
= + ч : 
| (Q, e, (178442 *845Q, e, *85,Q, e *(a4478,1855*84585,)Q, , 7 
б 271 Pe eo 172 2 
= 1 +а+ъ | 
(4.830) 
E aata d л 
f,, = UT (1.вз3за) 
letia Tt D 
But F must be a symmetric matrix, which 
implies that fo = Lor so that Equations (4.83b and c) imply 
that de,e, = e, making Q e a symmetric matrix, 1115 15 


ШОО сеп with the assumption that cross correlations are 


Porstant for all i, since when Q contains Rag i» Q 





Зе а 
] == ae = $ t е Е 
contains Re (4) Rag 1 Вов‘): Then the entries of 
become 
-а „- + 
( (158417855 а 2222214 е да бе e? 
TA Та | 
пи 2 O 
11 ] + а 4 6 
(340, „ "С -ајјаррћај рај ОФ е 8129 е 7 
112 5 15] 7 | 
нае (4.846) 
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| к | 
p. e, (1-84,*855 1187281282110 , ! 


mr a tD 





22 (4,846) 


Ehe entries of the vector G (or G') which appear in (4.48) 


can be written as 


(4.85) 


ee Qa. 


Since the object is to apply (4.84) ana 
(4,85) to the canonic realizations of Chapter III, it can 
be assumed the Qo gm О since for all realizations either 
e 


e(n) or e(n) is zero. Then 
gi 7 D o oU (4, дба) 


¿/(ltatb) (4,86b) 


рита 9 
2 21 e1 


When Equations (4.84) ana (4.86) are applied 


to the canonic realizations, then for the forms 


и = ше 1.872 
SM,TM: Te S А (4.872) 
x = Qe е /(1+a+b) ae 
1 
- 4.6 
Еро = 0 ( 7c) 
кг = =9@ _/(itatb) (4.674) 
a: 2 е Е 


RI 





BN TN: 11 


ша 


22 


61 


52 


since e,(n) is zero 


SM” TM”: 


EN ,IN : 11 


ke 


Да 


Че е (4. 88а) 
Ме a (4.885) 
0 (4.88c) 
(4 „вва) 
- Тен (4, 88е) 
for the untransposed forms, For the forms 
a + 2 рға жазы) (4,89а) 
Ee 2 
(-bQ +(1-b)Q +Q )/(1+a+b)(4.89b) 
т CI- o 
(200, e HA ии (4.89с) 
RES ғ. 0 (4.89а) 
pocos + ИРЕ (4.90а) 
Сазан, Жа в "Заре АА 
(4.90b) 


(-2(1*a4b)Q, , *(3*a-b)Q, , )/(1*atb) 
Бра eses 


(59007 





in 5^ ^ (4.908) 


These expressions can now be substituted 


Mito (4.48). This yields for 


lete) Q 2(cte)Q 
ES 8 2 
SM,SN,TM,TN: D E Cis a с (4,91а) 
Ni (l+a+b) (l+ta+b) 
Q + 29 + © 
е.е = |= Bere 
SM* , TM" : EE — LLL (4.91b) 
Ы (l+a+b) 
Q 
е е 
SN* tn’: DM (4.91с) 
Т (1+а+ь) 


Encre. Гог the T forms, c and e are replaced by c' and е! 
respectively. These may be written in another form by 


recalling that 


Q ser = Ее, (п-1)е (1) 1 = Ее. (п)е (1) ) (4.92) 


mor all i. Then for 


(ete)e, (n) 


2 
.— SM,SN,TH,TN: of = Етра + en)}“) (4.93a) 
t ont 2 e,(n) + esla) > 05 
SM',TH': E a Са 
2 
SN тие: Яс = E[(e5(n)/(1tatb)) ] (4.930) 
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When the filter is in a limit cycle condition, 
it is hypothesized that a deterministic mean-squared error 
analysis would suggest replacing Q e, by the mean-squared 
error (effective value). Since it is well known that the 
limit cycles of the quantization error sources are bounded 
by h/2 for round off, then for that quantization method the 
rms value of each source must be less than or equal to h/2. 
Thus the rms value of’ the composite errors еі» е, and € 


would have to be such that 





= 2 
e? 2 из В (4.9ШЦа) 
2 - % 
и М UK ан 
E дон (1.94) 
E 
Е. ne | Ц 
еје < ууу яг (4.944) 
Да 2 (oss 
Е IS vim . 
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where the bar indicates "mean" and H, and v are defined 
after equations (4.55) and (4.65a) respectively. Replacing 
the Q, e, in (4.91) by the right hand side of the appropriate 
equation from (4.94) results in hypothetical bounds on the 
rms value of the output error as given in Table 4-2. These 
are conjectural and a detailed analysis is suggested for 
EHCUPEe research. 

One example of limit cycles given by Jackson 
[10] has been examined using the formulae of Table 4-2. The 
rms value of the largest limit cycle in the example was 
4.24h. The formula in Table 4-2 yielded a maximum rms value 
of 4.44n. This example is offered not as proof but rather 


mee motivation. 


р. CONCLUSIONS 

Tables 4-1 and 4-2 list formulae for the output error 
variance for two extreme cases, uncorrelated error sources 
and maximum limit cycle rms values. It should be pointed 
out here that if any of the multipliers in the realization 
(not just in the transfer function) are such that it does 
not contribute quantization errors then these formulae need 
to be modified (by reduction of и ог YJ tO take this into 
account. For this reason, the expressions here cannot Бе 
directly compared to those previously presented in Тар1е 2-1. 
Another reason is that previous analyses have not included 
any dependence on zero locations. This dependence is expli- 
cit in the expressions for the untransposed configurations 


since operations involving zeroes tend to follow those of 


ea 





ARRAY 


TMATN 


SMoo »9Ngg » TMgTNo 


SN SN ТМ 


ЭМО 28M шары, 


102 


SM SN. SN ТМ 


B 20599929859 TM, 


02° 


SM, > ЗМ: Ny > SN 


SM4 4 ,SN44 





2 
(АУ) пах 





I ecke < h* 
1+а+р E 
г2(с+е)+3(1+а+) у2 h* 

1+a+b E 

г3(с+е)+2(1+а+р) -2 ne 

l+ta+tb СЦЕ 

1+а+р+е+е -Д 2 h? 

al 1+a+b ] 

-2(с+е)+0(1+а+ь) 92 h? 

l+tatb d 

2 


(се) #2 (1жать) 2 һ 
l+tatb 


d pos 
PA nn a n 
16 ғаз. aie 


1 2 pó 
dam? d 


TABLE !-2 MAXIMUM MEAN-SQUARED OUTPUT 
ERROR FOR HYPOTHESIZED 


LIMIT CYCLE ANALYSIS 
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poles in the flow between input and output in these realiza- 
tions, and thus are able to have a filtering influence on 
the internal noise generation. In the transpose configura- 
tions zero operations precede pole operations so the zero 
locations only influence the output error by the amount of , 
correlation involved between error sources. 

Àn interesting observation regarding the expressions 
for error bounds by the several investigators is the recur- 
rence of various factors. For example, the factor, A, in 


Table 4-1 can be written 


ЕРГО нъ) as) (1.95а) 
- (1-b)(1-a*b)(1*a*b) (4,950) 
= (Neyo OPS ИКЕ») (4,95c) 


The factor (1-b) appears in the bounds of Sandberg and 
Kaiser, Long and Trick, and Yakowitz and Parker and in the 
nimate of Jackson. The factor (1-|а| +5) appears in every 
bound of Table 2-1 for real poles and in the single bound 
of Jackson and also of Bonzanigo. The factor A appears as an 
intermediate result in Parker and Hess [52] and Aggarwal [48]. 

Equations (4.10) and (4.21) are the most general result 
of this chapter. They could be very useful as a tool for 
understanding the interrelation among pole-zero location, 


structure ana correlation in the accumulation of quantization 
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Exsors 1п the realizations of digital filters. Equations 
(4.48) are the general result for the realizations of Chapter 
ТТТ. It is suggested that similar results could be found 
for the ladder structures by formulating their algorithms 
via state space techniques. Another area of promise is in 
that of multidimensional filters. State space formulation 
atnese algorithms could possibly be expressed in terms of 
tensor analysis and the expressions of CHO) 
likewise. | 

Another area of application for these results is in the 
examination of realizations which perform quantization of 
signal variables after additions. In this case it is 
Suspected that accumulated roundoff error will be reduced 
since the maximum composite error at time n would be cre 
Quantization step for truncation or half a step for roundoff 
rather than the sum of the maximum error from several sources. 
Quantization could also be performed before multiplications 
allowing double precision signal variables to flow as far as 
Possible through the filter delays. There is often at least 
one path from the input to the output which would not require 
quantization. 

шистгезит Ев thus far still do not indicate which is the 
"best" filter configuration. The reason is that the nature 
of the correlation between error sources has not been completely 
examined. Extension of the work of Appendix B is an avenue of 


endeavor which would be useful to pursue. Once this has been 
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ccomplished, the "best" configuration for a given set of 
coefficients could be determined when the nature and 
statistics of the input signal are known. 

The results of this chapter can also be used in the 
determination of the optimal pairing of poles and zeroes in 


a cascade realization of higher order filters. 
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Figure 4-1b. Realization of Tag including composite error sources 
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Figure 4-2a. Non-Linearity of Rounding 


e (fane f()- [F1]. 






h 
2 


Figure 4-2b, Error Function for Rounding 
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V. DISCRETE FOURIER TRANSFORM FILTER 


A. INTRODUCTION 

In this chapter some relationships between the frequency 
Momain and time domain representations of signals in finite 
impulse response (FIR) filters will be examined. These 
relationships provide insight and motivation for developing 
anew convolution generator algorithm which yields the cir- 
cular convolution of a given sequence (Such as a finite 
impulse response) with sections of data taken from a continual 
data stream. Similarly, a new algorithm for performing the 
Discrete Fourier Transform of such sections of data is pre- 


sented. 
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mo) DEFINITION OF TERMS 
Recall from section IIA4 that the circular convolution 
of two sequences of length N was defined by the sequence 


B. Пп - 0,1,...,N-1 such that 


N-1 2 N-1 
= È h = È 
Yn meo м «n-m»? o mm (5.1) 
where {ui}, (а) n = 0,1,...,N-1 are the two sequences to be 


convolved. The notation <.> indicates modulo N arithmetic.! 
Subscripts (or superscripts) are used as spatial or positional 
Еп ои ега, whereas indices in parentheses will indicate 
time of occurrence or functional dependence. 

It was also shown that the same result (equation (5.1)) 
could be achieved via the equivalent operation of multiplica- 


tion in the frequency domain by defining: 


N-1 
DE © wwe = prrf{u.}] k = 0,1,...,N-1 52) 
k п-0 1 n 
N-1 ok 
ШО fh Ww = РВ РО k = 0,1,...,N-1 cor) 
n=0 
ШОН, Е бо (5.4) 





Lro find <x> = x modulo y, write x = r + py such that 
p = [x/y] where [-] indicates the largest integer less than 
Ше саца1 со (.). Then 0 < г < у and 
ШІ? - x modulo y = r = X - Ру. 


M 





Be 
Ук 7 OIDPTDU] no- 0,1,...,N-1 — (5.5) 
where W = exp[-j27/N]; (j = V) (ЕЕ) 


Note that wot is che principal wen РОО. 


The subscripts, k, above indicate the harmonic of 


interest, where the fundamental frequency is 
2 = 2т/МТ ПЕН) 


End T is the sample period in the time domain. 
me acion (5.5), if n is allowed to be outside the 
саси [0,N-1] then а Е 67 i.e., the infinite sequence 


ly, is periodic in n with period N. 


.- 00 
ШИЕ THE DIT FILTER 

ime сегт DET FILTER, as used here, refers to a form of 
a finite impulse response (FIR) filter presented by Gold and 
Jordan [53]. This form will be examined in detail to gain 
knowledge about the nature of the output, thus providing a 
stepping stone to the development of a new convolution 
meneretor algorithm. 

1. Formulation 

Consider the impulse response, h(n), of an FIR 
filter, which is zero outside the interval 0 < n < N-1. This 
function has a Z-transform given by 
N-1 


H(z) = E h(n)z (Seon 
n=0 
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which is the transfer function of the FIR filter. Define the 
sequence, е) и = О uch hat A, - h(n).* Then 
h(n) may be replaced by the Inverse Discrete Fourier Transform 


EN ts Discrete Fourier Transform, HH, as defined earlier. 


Then 
N-1 , N-1 N-1 Nel 
H(z) 2 I (y za? -l zr quy r q*z1)] (5.9) 
n=0 k=0 . k70 п-0 
N-1, , N 
Using the identity Ех = L—— — , yields 
п=0 шы 


12-0 М1 Н 


Н(2) = ше s loa Ul (5.10) 


since y "EN = exp[j?ık] m 

Equation (5.10) expresses the transfer function of the 
FIR filter in terms of the DFT coefficients, (HJ, of the 
Sequence Ша) = {h(n)}. This is the form utilized by Gold 
and Jordan to show that "finite impulse response" and "non- 
recursive" are not synonymous since, as shown in Figure 5-1, 
each term of the summation in (5.10) is a recursive filter; 
yet H(z) represents a finite impulse response filter. 


But there is more involved in this form than is 


immediately obvious. 


“The dual notation 15 used in an effort to maintain 
consistency in the subscript/parenthesis notation. {h(n)} 
and {ih} play a dual role here since (n, Jj indicates coefficient 


Mesition in the filter, while, as an impulse response sequence, 
(h(n)) is temporal in nature, yet they are identical in an 
ШЕН filter. 
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consider the terms, H,, in (5.10). Not only 


lm 
th з 
is Hy the k в пос of the DFT of {htn)} but it is alec 
E varue of the transfer function, H(z), evaluated at the 
ща th 


К М -root of 1, i.e., from (5.3) and (5.9) it can be 


seen that 


2. (5.11) 


Second, note that, since wk is, in general, a 
complex coefficient, the terms of the summation in (5.10) 
represent first order complex filters. Appendix shows 
that every first order complex-coefficient filter is also a 
song order real-coefficient filter. Furthermore, since 
ІНЕ eee poles of these filters lie on the unit circle 
Bethe Z—plane; that is, each term represents a digital 
resonator with frequency 
qu = KR = p (5.12) 
(Cf. equation (5.7)). In other words, the form of the section 
of the filter in 5-1 between point (4) and the cross-section 
- QI that of an oscillator bank and lew is a comb 
filter. The algorithm represented by equation (5.10) will 
henceforth be referred to as the DFT FILTER. The motivation 


for this phrase will become more obvious in succeeding sections. 
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ши Незропзе to a Finite Duration Input 

It is well known that the response of any filter to 
an arbitrary signal is the linear convolution of that signal 
with the impulse response of the filter. The relationship 
between the response of an FIR filter to a finite duration 
signal (of length N) has an interesting connection to circular 
convolution. A more complete connection is made if one exam- 
ines the response of the FIR filter to a signal comprised 
of two successive copies of this signal of length N (total 
length 2N). Then a method for finding the circular convolu- 
mon can be devised. 

The response of any filter, H(z), to any input signal, 


u(n) is given by the linear convolution: 


im 8g 


y(n) = 
m 


u(m)h(n-m) (5.13) 
When H(z) represents an FIR filter with an impulse 

response of length N, (h(n); n=0,1,...,N-1), and u(n) is zero 
Sept for n = 0,1,...,N-1, then by dropping all terms in 


(5.13) known to be zero (see Figure 5-2), 


п 
» u(m)h(n-m) 0«n « N-1 
m=0 
zu 
B. (5.11) 
у(п) = 5 u(m)h(n-m) N-1.< 1 < eN 
m=n-N+1 
0 otherwise 
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since only on the heavy black portions cf the abscissa in 
Figures 5-2 d and e are the products in (5.13) non-zero. The 
only time (5.14) yields a term of the circular convolution 
of the two sequences {h(m)} = ed and {u(m)} = (а) 
m= 0,1,...,N-1, is when n = N-1, since this is the only time 
when there are N products of terms taken from (u(m)) and (h(m)) 
in the summation. Then 
N-1 
y(N-1) = % u(m)h(N-1-m) (Ба Ша) 
m-0 

вас this is a circular convolution term is assured by the 
meet that, in (5.14), because of the limits on т, (п-т) = <n-m>. 
Buc this is only one term of the circular convolution. 
ШО па the other terms, it is necessary to resort to another 
technique. By passing two copies of the data through the FIR 
Entrer all terms of the circular convolution are obtained. 

Since circular convolution can be viewed as linear 
convolution if one sequence is replaced by its periodic 
Mitension, consider the response of the FIR filter to ап 
input comprised of two successive copies of the sequence 


mom) }, m = 0,1,...,N-1, i.e. 


u(n) O < n < N-1 
u(n) = u(n-N) N < n « 2N-1 (5.195) 
0 otherwise. 





The response to u(n) would be (see Figure 5-3): 


n 

Enemy mI) оо е пеш 
m=0 жа 

n 

5 u(m)h(n-m) мо 
m-n-N-*1 др 

y(n) = (Ба) 

2Х-1 

E u(m)h(n-m) 2N < n < 3N-2 
m=n-2N+1 Е 

0 otherwise. 


But in the second summation the limits on m imply 
that O < n-m « N-1 so it can be written immediately that 
n 


К 5 u(m)h(<n-m>)  N « n < 2N-1 
m-n-N-t1 


^ 

wn 
H| 
CO 

~ 


ши шеттмоте, because of the constraints onin mene 
Шомег limit on m is always less than or equal to N, so (5.18) 


бап бе broken into two parts, yielding 


N-1 n 
y(n) = x u(m)h(«n-m») + X u(m)h(«n-m») (5:19) 
m=n-Nt+1 m=N 


N En 5 21-1 


But since u(m) = u(m-N) for N « m « 2N-1 (from equation 


(5.16)), it follows that 





N-1 n 
2 u(m)h(«n-m») + E u(m-N)h(<n=m>) 
m=n-N+1 m=N 


y(n) 


N-1 n-N 
= E u(m)h(<n-m>) + E u(m)h(<n-m-N>) 
т=п-М+1 т= 0 
5-1 
= 2 u(m)h(<n-m>) N <n < 2N-1 (5. 20) 
m=0 
Since <n-m-N> = <n-m>. The N outputs, y(n), occurring for 
N «n < 2N-1 are exactly the N members of the set of terms 
resulting from the circular convolution of the sequences 
{u(m)} and {h(m)} m = 0,1,...,N-1. Thus the fact that circu- 
ват convolution can be viewed as a linear convolution if one 
Sequence is replaced by its periodic extension has been 
ehown írom a new perspective involving the FIR filter. 
But when one is processing continually received 
Mata, it is not efficient to stop and form periodic exten- 
mons Of blocks of data to achieve circular convolutions on 
пас block. 
Instead, the Fast Fourier Transform (FFT) has been 
used to perform high-speed circular convolution on blocks 
of data [ 15 - 19]. When these blocks are non-overlapping, 
MES is a highly efficient method and there is little need 
to overlap when the signal being examined is highly stationary 
Mia statistical sense. However, few real processes are 
highly stationary. To overcome this, overlapping blocks of 
data are processed. The question then arises, "How much 


Overlap does one need and what will it cost in computation time?" 
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The following development provides an algorithm which yields 
the results in 2N-1 operations for complete overlap of data. 
This development for obtaining the circular convolution of 
blocks of data with the impulse response of an FIR filter is 
Suggested by the symmetry property of convolution and the 
concept of periodic extensions of sequences. In equation 
(5.16) the periodic extension of {u(m)} was defined in order 
IA the circular convolution giweneby (5.20). -Why not 
use the periodic extension of {h(m)}? Then the circular 
convolution of {h(m)} with {u(m)} could be found by the 
Jinear convolution of one copy of {u(m)} with the periodic 
extension of {h(m)}. Some very straight forward methods of 
Being this are possible. But with additional insight, a 
metnod of accomplishing this ch the fewest number of summers 
and a simpler topology is possible. This is the object of 
the next two sections. The straightforward methods will be 
discussed after this development for comparison. 
3. ootate Analysis 

The previous section states the objective of finding 
an algorithm which provides the circular convolution of a 
sequence {h(m)} with one block of data from a signal u(n). 
im further suggests that some modification of an FIR filter 
WOuld be suitable to this purpose. However, in that section 
only finite length signals were considered. From this point 
EU the objective is generalized to that of finding the cir- 
cular convolution of {h(m)} with every block of data of 


length N taken from an arbitrary continual data stream, 
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mn). Examining the states? of the DFT FILTER will provide 
valuable insight into the nature of the signal quantities 
involved. 

Referring to Figure 5-1, define the states of the 
oscillator section as x,(n), k = 0,1,...,N-1 and let x(n) 


be the signal quantity at point (A) 3 помен 
хп) = u(n) - u(n-N) | (Sen) 


where u(n) is the input signal, which is assumed to be zero 
morn < о. 

For bookkeeping purposes, define the sequence 
{ч (а)), p = 0,1,...,N-1, (where n is now a time index), 


such that the elements of the sequence are 
u fn) = u(n+p-N+1) р = 0,1,...,N-1 (5.22) 


These sequence elements may be viewed as samples 
taken from a set of functions up (n) (see Figure 5-la for 
N= 4) also defined by equation (5.22). For example consider 
the sequence, {ч (п!)} = (3, 2, 1, 2) shown in Figure 5-85 
plotted versus p. The values 3, 2, 1 and 2 are the func- 


tional values of и (а), pon). u,(n) and uz(n) taken at 


1 
біте п! іп Figure 5-4a. Notice that Figure 5-4b is identical 


States are taken to be the signal quantites at the inputs 
to delay elements of the resonators of Figure 5-1. 
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mne region of the graph of u(n) (the top graph of Figure 
5-4a) indicated by a heavy black line on the abscissa. This 
heavy black line is the "window" for viewing the block of 
data referred to as tu,(n!)}. The time index n' marks the 
right hand edge of the window. The left hand edge of the 
window is then n'-Nt1. 

Next, define the sequence (8, (1) 1, к = 0,l,...,N-1 
as the DFT of Са (п). Note that the time index, n, associates 
these two sequences in that ts, (n)) is the DFT of the block 
of data in the window at time n, which is by definition 
(и Сп) 3.3 Specifically, 

N 


з (п) = E 
k р-0 


T N 


1 
pk 

u (n)W = 

p ) 


Y u(ntp-N41)wPE ERO 
p=0 
The relationship between the states of the DFT FILTER, 


x(n), and the sequence {в (п) } can now be found. 


From Figure 5-1, it is seen that 
x, (n) - W Fx, (n-1) + x(n) (5.24) 


By applying (5.24) iteratively (m-1) times it is 


Seen that 


m=i 
x, (n) = Мх (п-т) + E ИРЕ, (n-p) 05825) 
р=0 





alt can also be seen that {u_(n)} is exactly the set of 
ис јез 5богед in the N-delayPof Figure 5-1 at time n. 
END is the quantity residing closest to but not yet arrived 

at the summer. The quantity at the summer is u(n-N) = ч0(0-1). 
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Assuming zero initial conditions, i.e. x, (n) - 0 
Гог п < -1, then for m - ntl, (5.25) becomes 
b ns 
x(n) = Z W** x(n-p) (5.26) 


Substituting (5.21) into (5.26) and assuming n > N-1, 


2 -pk i -pk 
х (п) = Z u(n-p) W P@ _ у u(n-p-N) КР (5.27а) 
р-0 р-0 
N-1 n n 
ТЕ) Wee + x u(n-p) w PES еи 
р+0 pzN p=0 


Performing a change of variables in each of the three summa- 
tions of (5.27a), i.e. q = Nep-l, q = p-N and q = p respectively, 
N-1 n-N 


x, (n) = £ u(nt+q-N+1) y (9-N+1)k + X u(n-q-N) М 
ал0 а=0 


-(q+N)k 


u(n-q-N) np cs (БОО үр) 


0 


| 
Шке 


9 


Combining the second and third summations of (5.27b), realizing 
Me y CI FWXK 2 yk violas 


n 
-pk 
x, (n) - "s X u(ntp-N*1) М E » u(n-p-N) W p 
p=0 p=n-N+1 
(52270) 


The summation in the first term of (5.27c) is, by 
№.23), s, (n). The second summation is zero by the assumption 


Beat ulm) = (С form < 0. So 
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ПЖ we s, (n) (5.28) 


Equation (5.28) shows that for any n > N-1, the y un 


state of the DFT FILTER is proportional to the К harmonic 
of the DFT of tne sequence, tu,(n)}, that is, of the sequence 
consisting of the present input and the past N-1 inputs. (In 
ШОО, this is true for n « N-1 if one considers Sequences 
with leading zeroes.) Thus (5.28) provides the motivation 
behind the name DFT FILTER for the algorithm of equation (5.10) 
mo Figure 5-1. 

The state response of the DFT FILTER for two specific 
types of input sequence will now be considered as examples 
Бо illustrate the operation of the comb filter section. 

a. Example 1: Unit Pulse Response 

bet win) consist only of the Unit Раш елитној 
This pulse is sent directly to the resonators as x(0) and 
Шатев their oscillations. It is also stored in the N-delay 
of Figure 5-1. (Figure 5-5 shows the result of this input 
for N = 4.) These oscillations are Such that x, (n) = a 
eo 1,...,N-1. 

Brom (5.22), ug (1-1) = u(0) = 1 and Це - 0 
mer p = 1,2,...,N-1. So fu, (N-1)} is an impulse sequence 
with the impulse at p = 0. On the other hand, xy (N-1) 15 
Bauval to ут = E as expected from (5.28) since the DFT 


k 
of an impulse sequence is {s,} = (1). The factor W can be 
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viewed as a phase delay of ?ık/N since ү“ - exp[-j2Tk/N]. 
This phase delay is evident in Figure 5-5. 

Recall that the original pulse was stored in the 
N-delay and now impinges on the resonators as a negative 
pulse at time n = N since x(N) = u(N) - u(0) = -1, But 
х (№) = ws x, (N71) + x(N) = 0 and so for n > N the states 
ШІ сће DFT filter will be zero. 

In effect then a unit pulse input causes the 
resonators of the oscillator section to oscillate for one 
period, n = 0,1,...,N-1, after which time they are turned 
off. Since an arbitrary input is simply a string of weighted 
and delayed impulses and since linearity applies, the state 
response of the DFT FILTER to an arbitrary input is the 
Superposition of a string of weighted and delayed oscillation 
periods. 

Before proceeding with example 2 some additional 
insight into the intimate relationship between the states 
Se the DFT FILTER and the DFT of sequences will be helpful. 
Consider equation (5.23) which by a change in variable 
(p is replaced by p-1) can be rewritten as 
N 


А (р-1)к 
s, (n) = ра М era) 


БЕН (5 22) іс follows, in general, that 


ц. (п) = u((n-m) * (p*tm)-N*1) = ш, (п-т) (5.30) 


р 


are 





Mor any m. So (5.29) becomes 


К N uy (n-1) РК 
p 
N 
= ( 


p 


5, (n) 


-1 
= М È u_(n-1) РК + u,(n-1) wNk -= u ШЕПТІ! 
=) Р М 0 


МТ (s,Q-1) # ша (а) - ша (а-1)) 


TS ts, (n-1) а = оспа) 


МТ {5 (п-1) + х(а)) (5.31) 
Bubstituting (5.31) into (5.28) yields 


( Е К 

x, (n) s, (n Пр ар раса) 5532) 
Bauacion (5.32) says that x, (n) is also equal to 

the DFT of tu, (n-1)) plus an additional term. By applying 


(5.31) iteratively (m-1) times it can be shown that 


т- 1 
в (п) = ММК 5 (п-т) + z WTCT Cng) (Б р 
Р к 


Substituting (5.33) into (5.28) yields 


m-l 
үл С (5.34) 
qz 


-(m-1)k 


= W = Үр + 
x, (n) W s, (n m) 
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So the states of the DFT filter at time n contain a phase 
shifted version of the DFT of previous sets of window data 
pius some additional terms. Example 2 will demonstrate 
Pots With a periodic input signal. 

b. Example 2: Response to Periodic Input 

Let u(n) be zero for n « 0 and periodic with 

period N for n > 0 (See Figure 5-6). Then x(n) (point (A) 
in Figure 5-1) consists of only one period of u(n) while 
the window of data stored in the N-delay after n = М-1 
mooks like a continually progressing circular shift of that 
period. The DFT's of such circularly shifted sequences 
is known to be identical except for a phase shift. Since 
x(n) = 0 for n > N the resonators will continue to oscillate 
Eun cns numbers of constant magnitude but changing phase 
at the rate, 2"k/NT radians/seccnd. This ls predicted by 
equation (5.34) when m is not so large that x(n-q) is a 
term from the single period of data which enters the 
што | тасот рок топ of the DFT FILTER. 

In the next section significant use will be made of 
Bue fact that x, (n) contains some information about the 
DFT's of any number of sequences preceeding fu, (n)} in time. 
maat this information is contaminated by phase shifts and 
additional terms will not prove insurmountable. 

ШОО ри Analysis 

The relationship between the states of the DFT FILTER 

ana the DFT of the input sequences provides, as shown below, 


the knowledge which makes possible the extraction of complece 


214 





information about the circular convolution of every block 
of data of length N with the impulse response of the FIR 
miter. 
According to Figure 5-1, 
N 


jeu 
y(n) + N 2, 
k=0 


1. 
нү x, (n) (5.35) 


But in the previous section (equation (5.28) it was 


Shown that x, (n) was proportional to s, (n). Thus, 


ш 
= L = 2. | -(N-1)k 
y (n) N > Hy W S (n) ХН 5 (а) 


which is, from equations (5.5) and (5.1) the (N-1)th term 
of the circular convolution of the sequence th, } with the 
sequence Си (а), since ts, (n)) is the DFT of und}. 


So the output may be written as 


N-1 N-1 
у(п) = 2 и мл» 2) h(m = Z ua) h(«N-1-m») (5.37a,b) 
m=0 m=0 
or 
N-1 N-1 
y(n) = 2 Uy y (m0 him) = 2 u,(n) h(N-1-m) (5.370,83) 
m=0 m= 0 


since «N-1-m» - N-1-m for the values of m allowed in the 


summations. These four expressions are equivalent. 





However, if equation (5.32) is substituted into (5.35), 
пеп 


N 1 №-1 
у(п) = 7 ОН 5, (п-1) АН 


К х(п) (5:59) 
k=0 k=0 


be first summation is the zero term of the circular convolu- 
tion of (n, ) with fu, (n-1)}. Since x(n) does not depend on 
k, the second term of (5.38) is merely x(n).h(0) since 

сін ОСУ cetinitvion, the zeroeth term or tne Грешно 


(Hj. Therefore, not only is y(n) the (N-1)th term of the 


(и (а) ) convolution, it is also possible to extract the 
zeroeth term of the convolution with fu, (n-1)} by subtracting 
EE h(0) from y(n). In fact, N-1 more convolution terms 


can ve extracted from y(n), as will be shown next. 


emeetitutine (5.34) into (5.35) yields 


N-1 N-1  m-1 | 
k=0 k=0 а=0 
N-1 | mel 
al zm s (nm) WIE 4 (mg) nla) (5.39) 
k=0 q=0 


the first term being the (m-1)th term of the (а (а-т)) 
convolution. For notational purposes define Ya nK) as Che 
mth term of the circular convolution of the sequence 

tup (Kk) ) with thy}, where the index n indicates when that 
term is available from y(n). With this notation then, 


equation (5.39) says 
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У 16050-1) = y(n) - "P x(n-g) h(q) (5 105 
a=0 

This algorithm is shown in Figure 5-7, using the DFT 
FILTER form of H(z). Figure 5-8 uses the FIR canonic form 
of H(z) and utilizes the normal delays of that form in 
conjunction with a direct feed-forward path to form 
x(n) = u(n) - u(n-N). 

The algorithm of equation (5.40) might be called a 
"circular convolution generator." It is most efficient when 
one must calculate all N convolution terms for every subse- 
quence of length N of the “infinite" sequence (u(n))g. In 
ас case, it requires, on the average, roughly two multi- 
plications per term (i.e., 2N multiplications per convolution) 
instead of the N multiplications per term required to perform 
each convolution separately or the (log, №) + 1 multiplice- 
tions per term required to accomplish each convolution via 
Ene FFT. 

It may be noted here that it is not necessary to 
 пегасе (һе term yy. .3 (n 5n-N) since this term is available 


as y (n-N) (n-N;n-N). This saves one delay, one 


a 
 ИБірізсасіоп апа two additions but is insignificant for 
large N. 

To illustrate when the terms of the convolution of 
a particular sequence (u, (m)} become available at the 


MiGpucs of the circular convolution generator of Figure 5-1 


or 5-8, Figure 5-9 shows the time base of the N outputs. 
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The points connected by the diagonal line represent the N 
terms of the circular convolution of (и (т)}. The last 
output is the redundant term Ум Q5n-N). 

Section B alluded to some straight-forward methods 
of performing this same operation. Figures 5-10 and 5-11 
are two such methods. The sequence tu, (m) } is stored in 
the first N delays of these algorithms at time n = m, there- 
fore these delay elements are equivalent to the N-delay 
Капас positions of the DFT FILTER. So by the notation 


th term of the circular 


convention, y, (M3k) represents the m 
convolution of tne data in the first N delays at time k and 
available at time n. 

Figure 5-10 requires 2N-2 delays, 2N-1 multipliers, 
and N(N-1) summers not including the redundant operations 
Bt the right end. The worst thing about Figure 5-10 is 
EE Chorrifying number of circuit crossings required. 

Figure 5-11 reduces the number of summers and circuit 
Serossings by recognizing that to achieve each successive 
EnCDput term moving from left to right it is only necessary 
bo add one product and subtract another. Thus Figure 5-11 
requires 2N-2 delays, 2N-1 multipliers and 3N-5 summers. 
Again the redundant operations have been removed from con- 
Exderation. But there are still N“ = ell eae Glance 
EPOsSsings! 

Figure 5-12 is a rearrangement of Figure 5-80 which 


Ras no Merc crossings.” This arrangement of the circular 


convolution generator requires 2N-2 delays, 2N-1 multipliers 





ч 


2This is particularly adventageous in Lol technology . 
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and 2N-2 summers. v is believed that this is the least 
number of operations required to obtain the objective. 
Obviously the topological circuit crossings cannot be 
reduced further. 

Шив шаг the portion of Figure 5-12 to the left of 


cross-section Q - (с) гоша be replaced by a DFT FILTER 


but to no apparent advantage at this time. 


Pee MODIFIED DFT FILTER 

In Section C3, it was shown that the states of the DFT 
Ж ШЕТЕН were proportional to the DFT of the input dence 
tup (n) j. It is also possible to design the oscillator 
section so that the states are exactly the DFT of tu, (mn), 
mars will lead to a minor refinement of the convolution 
generator. It will also lead to some observations that will 
be pertinent in Section E where a new DFT algorithm is 
derived., 


Taking the Z-transform of equations (5.24) and (5.28), 


Х (7) _ К 
НЕ EW Se) ке ом 
k 1. wok „71 k 
И) 
where 
X, (2) = Zix, (n)J (5.41а) 
X(z) = Ge (5.410) 
S (z) = Z(s, (n)j (5.410) 
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The right-hand equality of equation (5.41) may be 
rearranged to provide a transfer function between X(z) and 


5,(2), і.е. 


5 (2) wok 
X(z) 7 тик (5.42) 


The right hand side of (5.42) represents a digital resonator 
with a feed-forward gain of WR and a locpmezın of in | 
Figure 5-13 shows a modified DFT FILTER where (s, (n)) are 
we States of the oscillator portion. 


Designating the output of the modified DFT FILTER as 


Ши), the transfer function of this form is 


DININ E 


siZ r „ К (5.43) 
k=0 1 - W 7 


In the form of an FIR transfer function, (5.43) becomes 


N-1 
G(z) = 1 h(n+1)z" (5.44) 
n=0 
where h(N) is by convention equal to h(0) to be consistent 
with the concept of periodic extensions of sequences. It 
would also be possible to replace h(n+1) by h(<n+1>) which 


is consistent with the "modulo N" nature of circular convolution. 





Omnis form is related to the Goertzel algorithm [ 4] for 
finding the DFT (see Section IIAÀ ). 
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Upon examination of the output of this [ident ES 
Been that 


y N-1 
у) ет ЕН, эк) (5.45) 


mech is the zeroeth term of the circular convolution of 
(ы (а) } with th}. 


Ми са ји equation (5.33) into (5.45), 


N-1 N-1 m-1 | 
роп) = 5 У Н 5. (п-т) wom 4 а L Ho E u (ng) 
k=0 k=0 q=0 
m-1 
= у (n;n-m) * 2 x(n-q) h(q+1) (5.46) 
m o 


Mees Snows again that the output contains sufficient informa- 
Enn to allow the extraction of N convolution terms in a 
Similar manner as was done previously. In this case 
m-l 
y,(Q;n-m) » v(n) - E x(n-q) h(qt1) co.) 
q=0 
This algorithm is shown in Figure 5-14 with G(z) repre- 
sented in the FIR form. G(z) may also be implemented, of 
course, in the modified DFT FILTER form. In either case, 
Figure 5-14 may be referred to as the "modified circular 
Envolution generator." 
The modified convolution generator has no particular 
advantage over the one derived in Section C except that the 


temporal order in which the terms of a convolution associated 


22] 








with a particular sequence (и (а) become available is also 
the ascending order of the index, m, in (y (n;k)) , 
m= 0,1,...,N-1. Furthermore, this result occurs without 
relying on the redundant multiplication noted in Section C. 
Note also that instead of the apparent circular shift in the 
magex, m, there is, in Figure 5-6, a circular shift of the 
mex On the coefficients, п, and that the redundant multi- 
plication would be by a second ho which does not appear. 
Extending this idea, it is then possible to extract the 
Benvolution terms in any circularly shifted order by a 

Mitable circular shift of the his. It will be shown in 
Section F that by relying on the symmetry property of con- 
volution it will be possible to extract the convolution term 

EN Uv circularly shifted order by a suitable circular shift 

of the pos. This will have a significant impact on a possible 


meplication of the circular convolution generator. 


Е. A NEW DFT ALGORITHM 

In the previous sections, the primary interest was in 
пе domain extraction of circular convolution, though 
frequency domain techniques were resorted to in order to 
man insight into the properties of the available signal 
Muantities. In this section, a new DFT algorithm is derived 
by making use of some of those same insights. 

In Section D, a transfer function between X(z) and S, C2) 
Was derived which provided the spectrum of the sequence 


tu, (n)] Mamie output of the oscillator section (cross 
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section (8) - (в) ОТ Figure (5-13)). "This "part of the 


algorithm is not new. It has been derived previously by 
Dillard [54] as a means of generating spectra 

or spectral components of a series of sequential sequences 
бы (0) 2, as defined in Section C3. This algorithm may be 
described as a moving (sliding) window DFT or uncoupled DFT. 

As a moving window DFT, it may be observed that one may 
feed one sequence of length N into the resonators and obtain 
mac entire spectrum of that sequence in ме multiplications, 
the same number required to perform the DFT in its definitive 
Borm. However, the algorithm requires only N multiplications 
per spectrum, on the average, if one is interested in full 
Exerlap of the moving window. Further savings are realized 
if one is interested in only a few spectral components, since 
the algorithm provides each harmonic independently of all 
его, і.е., it is an uncoupled form. 

This algorithm has one major disadvantage which has not 
Бееп expounded or examined to date. Theoretically this 
algorithm provides the exact discrete spectrum of each 
sequence, (4,60), independent of every other (ы (а). Нои- 
ever, when finite precision arithmetic is used, significant 
errors destroy that idealistic situation. 

First of all, it is apparent from studies of possible 
pole locations in second-order, finite word-length, real coeí- 
ficient digital filters [ il ], that no method has been de- 
vised (or can be. devised) to place the poles of a set oi such 


E ters exactly on the unit circle in the Z-plane and evenly 


eS 





ИС ЗЕЕ frequency (except for the set іг|2- 1, 3, 21, -j}) 
because of the finite length of the coefficient representation. 
This means that there will be errors due to the fact that the 
transfer function realized is not the one required. From 
another point of view, the problem is tnat "gs in the defini- 
шоп ог the DFT is not representable, in general, in a finite 
шета length. So this difficulty applies to any realization 
moa DFT algorithm. 

ЕЕсепайу, іп а Гіпібсе word length algorithm; it ds 259 
mecessary to quantize the results of multiplication because 
of the extended word length that operation implies. This 
Em Pculty, too, applies to any realization of the DFT. 

But. the third and most significant problem to be faced 
Mere is the recursive nature of the comb filter section of 
Mesure 5-13. In effect, it magnifies the problem of quanti- 
шастоп агбег multiplication. 

Recursive digital filter sections have been extensively 
Studied in the literature [ 39, 51 ] and in Chapter IV of 
this work with respect to error generation and propagation. 
particular, recursive filters are highly subject to the 
generation of limit cycle conditions, especially when the 
Dole locations are close to the limits of stability. In 
mae case of the DFT FILTER the pole locations are exactly 
On the unit circle, the stability boundary. It should suisse 


here to say that recursion should be avoided whenever possible 


mI practicable. ` 
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There is a set of complex coefficients, however, which 
do not contribute errors of the type described above. These 
are the coefficients (1, 3, -1, -j). If a transfer function 
contains these and only these coefficients, there will be 
no shift in pole location between the required and the 
realized filters. Furthermore, there will be no error intro- 
duced by quantization after multiplication by these coeffi- 
cients and hence no recursive propagation Of- -such errors: 
Емеп іп a recursive filter. This seems to be a trivial idea, 
But is importance comes to light in view of an observation 
by Hess [11] regarding the realization of digital resonators 
ОГ а desired frequency. 

meso observed that, since the frequency of a resonator 
is dependent upon the sample period, any frequency can be 
generated. By choosing E equal to an integer multiple 
72, only multiplications by 1, j, -l or -J are required 
mo have an oscillator with frequency Ша as shown below. 

It was shown in Section D that the DFT of windowed data 
wna be achieved through an appropriate arrangement of 
digital resonators whose poles were located at the associated 
harmonic Ргеацепс ез. An algorithm will now be developed 
which exploits the relationship between frequency and sample 
period. 


Consider the complex resonztor with the transfer functron 





H (z) = ze: (5.18) 





H(z) has a pole at 


Ри = ехрГіт/21 (5.49) 


Equating the right hand side of (5.49) with the exponential 


form of z on the unit circle: 


A exp[jw¿T] (5.509 
yields the resonator frequency 


— ТІ d 


Furthermore, one may aiso observe that this poie is ас 


which means that this resonator would generate the DFT 


frequency 
ыссы (5.53) 


Ши this is exactly w, in equation (5.51). 


а 
Now, if T were to be replaced by T! = 2T then this pole 
location would correspond to half the original frequency. 


However, it is not part of the approach at this point to 


consider changing the sampling rate, since it is assumed 
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that data 15 being received at a fixed sampling rate. On 


the other hand, consider the complex resonator 


Но (в) = —= (5.51) 


l- jz 


and evaluate this function on the unit circle, Thus 


ed 
he AN 1 д 1 
Н (е Ша -ju 2T i je 17 (555% 
Je l- je 


This is exactly H, (2) evaluated on the unit circle when Т 
is replaced by T' = 2T. The sampling rate had not been 
Smanged in equation (5.55) but the effect is similar. Іп 
mact, if H, (2) were used in a desampling mode (where only 
every other output is considered) the effect would be the 
same. 

But now consider the poles of H5(z) which occur at 


1 2T 
7 = EN 2 = +[031/232 = Е "EE (5,56) 


Mie could achieve the effect of a single pole at 


N 
1 


EE by modifying H5(z) to include a zero at 


ехрід5т/4|. Thus 


N 
| 


j57/4 _-1 -5N/8 _-1 
ше = = Loto z (5.57) 


Шанг“ I 


al 





H (z) = 
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The pole, 7 - ехрГіт/!1 - wN/8 ls one required in any 


DFT FILTER algorithm which has N equal to an integer 
multiple of eight. 

Immediately obtainable from H, (z) is another pole required 
ENNMultipie-of-elght algorithms; i.e., z = и-2М/8, Ши fact 
all eight frequencies for N = 8 are obtainable from the set 
ши Сгайзтег functions: 


-5N/8 -1 
jJ Ri ^ 
ul) = a = 0100) = He) , (5.588) 
EP. I 


ЕА c 
Соцув( 2) - т же = G, (z) (5.58c) 


-2N/8 -1 
44 -W Z d 


JM 2) 


ы wo IN/8 гБ 8 
Gan /g 62) m poc 6-02) (Ss е) 


м-3М/8 ,-1 


1 - Z 
G Ca) А G- (Z) (5.587) 


M S 


2 = 
Gon 7g (2) IL an Ga Cz) (5.58g) 


wo -1 
= "= = в, (=) (5 Бош 


Симув( 2) EE 


Where the subscript indicates the harmonic number of the 
associated frequency. That is, the subscript corresponds 


БОК іп Figure 5-13. 





Figure 5-16 depicts an algorithm which makes use of this 
development for М = 8. As in Figure 5-13 the input must 
first pass through the section (1 - P before reaching 
bhe resonators. Also, recalling that each resonator in 
Figure 5-13 had a feed forward gain of uA each forward 
path in Figure 5-16 includes this scale factor which was 
Hot included in equations (5.58a-f). 

Figure 5-16 is a rather complicated algorithm for an 
№ = 8 DFT. It also appears to require twice as many 
multiplications as the algorithm of Figure 5-13. Some 
further refinements would make this a more desirable form. 
Examination of equations (5.58) yields these refinements. 

Consider equation (5.58h). Multiplication by ТЫ 15 
puvial since we = 1. In addition, the denominator of 
M5 560) is factorable into (1 + gU (1 - 271) which means 


(5.58h) can be rewritten as 





Ji 
О E Door) 
4 DE да 
fe (5,50¢) a similar result occurs. Here wo ИВ = =l, 
со that 
t _ 
С (2) = ae „71 (5.590) 
жш (652552). ноиб - j; thus 
1 - dz 1 (5.59с) 
G = = non = 5 С 
Ще 2T 
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-2N/8 . 


ЕЮКСГІП (5.5050), М - 1 пат 


il 


Qr iz) = == 
6 1 - 3271 


(5.594) 

Equations (5.59a-d) denote resonators with a single delay 
element, rather than the two delay elements in equations 
(5.58 c,d,g and h). The cancellations which were used to 
derive equations (5.59a-d) from (5.58 c,d,g and h) are not 
possible in equations (5.58 a,b,e and f) without reintro- 
ducing non-trivial multiplications into the recursion loop. 
(A non-trivial multiplication is one where the coefficient 
is other than 1, j, -1 or -j.) 

Another simplification involves flow diagram manipulations 
to be performed on the zeros associated with the resonators 
of equations (5.58 a,b,e and f) including the scale factor 


-kK 


W associated with the appropriate harmonic. 


Figure 5-17a shows the zero paths associated with 
equations (5.58 a and b) after observing that w- 28/8 is 
-N/8 


equal to -W (In fact, this is tEUCOUndepemwent о нв 


porue of N, since 


и“ 748 « euo (-322/N) (-IN/8)] = ехр[јт] = -1.) 
(5.60) 


It is now possible to remove the multiplications by 


и-К/8 which are closest to the output in Figure 5-17a by 


> -N 
Placing a multiplication by W /8 in each input to the 





Summers aS shown in Figure 5-17b. Then one can consolidate 


Sensccutvive multiplications by w-N/8 to obtain a multipli- 


~2N/8 
cation by W / which is actually a multiplication by |] 


as shown in Figure 5-17c. The next manipulation is to 
remove identical multiplications which occur in all paths 
leaving a node by placing a single multiplication before 
that node. Figure 5-17d shows the final arrangement of a 
flow chart for equations (5.58 a and b) which is equivalent 
to that shown in Figure 5-16. (Residual multiplications by 
-1 have been moved to the right and made part of the summer 
(subtracter) operation.) A similar manipulation is possible 
for the zero paths of equations (5.58 e and f). 

Figure 5-18 depicts an equivalent flow diagram to Figure 
5-16, incorporating 211 the above-mentioned modifications. 
(Note, only non-trivial multiplications are drawn inside a 
Mare. All trivial multiplications, i.e., multiplications 
By 1, j -1 or -j, are shown inside circles.) Three major 
Meoperties of this figure should be observed. 

Horst, all non-trivial multiplications have been left 
En the form wkn/B because in this form they are independent 
of N. On the other hand, the subscripts on s, (n) аге also 
shown in this manner to indicate their dependence on N. The 
usefulness of this fact will be seen later. 

Second, Figure 5-9 15 а somewhat uncoupled госп. He 


four single-delay stages are completely uncoupled and tne 


two double-delay stages are pairwise uncoupled. 
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Third, there are no non-trivial multiplications in the 
recursive part of the algorithm. Therefore, no recursive 
error propagation will ensue. Errors generated in this 
operation are only in the output zero-paths and only effect 
the DFT of one window of data. 

A final observation on this specific example is that the 
number of non-trivial multiplications involved is reduced 
from what was required in Figure 4. There are six nontrivial 
multiplications required in Figure 5-13 for N = 8, while 
only two exist in Figure 5-18. A first estimate on the 
mmber of multiplications for the FFT is N log, N = (8) (3) = 
ey, but through manipulations of flow charts similar to that 
merrormed on the zero paths; the FFT algorithm can also be 
meauced to two non-trivial multiplications. So nothing nas 
been gained in terms of non-trivial multiplications over 
the FFT, but some degree of uncoupling has been achieved. 
This will shown to be an advantage in succeeding paragraphs. 

When N = 16, four resonators With four delays each) are 
required to provide the 16 poles necessary for the transfer 


шцис 1015, 


ит КК/ 16 


G Qr ата ШЕ (5.61) 
kN/16 i oe 


Again the feedback coefficients for the four resonators will 


ШЕ 1, Jj, -1 and -J. 
Ше. 15 песеззагу in this case to cancel three poles of 


each resonator by zeroes to provide the single pole of (5.00 


adz 





In the process of doing so 1% сап be observed that the 
resonators with feedback coefficients 1 and -1 reduce 
exactly to Figure 5-18 with N = 16. (Notice that N = 16 
does not change the coefficients but only changes the sub- 
Scripts on the outputs, since now W - exp[-j27/N] » exp[-j27/16].) 
This reduction occurs through the same reasoning by which 
Equations (5.59) were derived. 

The remaining two resonators and associated zero paths 
are shown in Figure 5-19. When the appropriate flow diagram 
manipulations are performed, Figure 5-19 reduces to Figure 
5-20. Figures 5-18 and 5-20 then, taken together form a 
complete DFT algorithm Гог М = 16. The outputs of Figure 
5-18 are the even harmonies and those of 5-20 the odd ones. 
ши Сре algorithm of these two figures taken together were 
бре built or programmed, there would be available to the 
user DFT algorithms for N = 1, 2, 4, 8 or 16. The only 
parameters to be chosen are: 

(1) the length of the N-delay at the front end 
(2) the number of resonator loops to take advantage of 
Gaye che subscripts of the outputs 
But all of these parameters depend on N and can be incor- 
porated in one switch or program statement which indicates 
the value of N. Notice also that it is only necessary to 


store two non-trivial complex coefficients since 


y-9N/16 ,. 4738/8 , j-N/* (5.622) 
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w-9N/16 . -N/16 


(5.625) 
апа 
-2М/16 -М/8 

1 - иу (5.630) 
Borthät only wN/8 = exp[ljt/4] ana y 3N/8 = exp[j3n/4] need 
stored. These can further be reduced to two real 
coefficients since 

-N/8. _ -3М/8 
ІМ ] = тм 1 (5.64) 


Another reduction of the algorithm is possible if only 
meal data is being processed. Then all the resonators with 
-) аз the feedback coefficient can be removed since the 
Outputs associated with them are available as the complex 
menjugates of the outputs of the resonator with the same 


number of delays and feedback coefficient, 4); e.g. 


524/16 0) = 55/16 m) (5.65) 


where * indicates complex conjugate. The output 520/16? 


ШЕ іп Figure 5-20b while Son/16 is in Figure 5-20a. 
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Algorithms for N equal to a power of 2 are possible 
for all higher powers as well. Figure 5-21 shows an eight- 
delay stage with feedback coefficient, j, required for N - 32 
and higher powers of two. Figure 5-22 shows the equivalent 
stage after reduction of the number of non-trivial multipli- 
cations. Figures 5-23 and 5-2! show the sixteen delay stage 
required for N = 64 and higher powers of two. 

Only the stages with feedback coefficient j are shown 
since the form for the stage with -j is the same except that 


-N | : 
/32 in Figure 5-22 are replaced by the same 


-N/61 


all powers of W 


powers of V and all powers of W in Figure 5-24 
are replaced by the same powers of и" 3М/64. 

Taking Figures 5-18, 5-20, 5-22, and 5-2! together as 
one algorithm, there are available the routines for finding 
ШЕ DFT with N = 1, 2, 4, 8, 16, 32, and 64, and the only 
parameter that needs to be changed is N. 

Another advantage of this algorithm is that the semi- 
uncoupled nature of it allows one to reduce the number of 
calculations required when only certain subsets of the har- 
monies are required. For example, Suppose one is interested 
ШІу іп the odd harmonics of real data for N = 64. Then 
only the algorithm of Figure 5-24 need be performed (inclu- 
ding ле feedforward о path of Figure 5-18). If only 
every other even harmonic beginning with k = 2 for N = 64 
is desired, then Figure 5-22 ylelds all the речи сео 
6 


р -64 : 
information, when (1-z ) precedes it. 





If the data being processed are not from a continual 
data stream, but rather are Just N pieces of data, then the 
delay = ls not necessary. The data can be fed sequentially 
to the resonators as long as the delay elements have been 
cleared (set to zero initial conditions) before entering 
the first datum. The calculations on the zero paths need 
only be performed at time n = N (if it is assumed that the 
first datum is entered at time n = 1). 

From the previous paragraph and examination of the 
resonators of Figures 5-18 through 5-24, it becomes clear 
that this method of finding the DFT is essentially different 
from other methods since the first operations performed on 
the data involve adding all data, every other datum, every 
Bourth datum, etc., after appropriate multiplication су 
mewers of 1, j, -1, -j. That is, if this method were performed 
in parallel from the start rather than with a sequential data 
feed, Figure 5-25 would result for N = 8. This also дешоп- 
Btrates that the sequential data feed is easier to implement 
in this algorithm; and,with the inclusion of the N-delay, the 
DFT algorithm of Figure 5-18 allows complete overlap of 
successive windows of data with only 7 additions required 


to set up the data rather than the 40 additions of Figure 


E-ob. 
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Е. CONCLUSIONS 

The convolution generator and modified convolution 
generator algoritnms produce the circular convolution of 
every block of data in a continual data stream with the 
(finite) impulse response Of an FIR filter. The algorithms 
are simpler in the number of operations to be performed and 
in topology than the definitive form of generating circular 
Eonvolutions. 


The men 


impulse response of the N outputs is & сору of 
the original FIR response circularly shifted m times and 
then delayed m sample times. Because of this the algorithm 
could be used to provide N matched filters in one algorithm, 
where the set of transmitted signals are the N circularly 
shifted versions of a sequence of length N. The maximum 
signal-to-noise ratio for a given transmitted signal will 
occur at the intersection of the diagonal line of figure 

БЕ 15 and the output line corresponding to the circular shift 
of the transmitted signal. 

The new DFT algorithm presented here allows the program- 
mime Or wiring of a DFT of length N = 2 to be used for 
performing the DFT of lengths a lower power ЈЕ Ewo рупе У 
changing the parameter of length. The algorithm is designed 
for sequential data feed which allows complete flexibility 
of overlap in processing sections of data from а сорта! 
data streem. The zero path operations need only be performed 


when the last datum of a particular section of date has 


entered the resonator loops. No quantization errors occur 
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in the loops providing complete cancellation of previous 


data not in the section under consideration. 
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Figure 5-2. Convolution of Two Sequences of Length N 
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gure 5-2. Convolution of a Length-N Sequence with the Single Period 
Extension of another Length-N Sequence 
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Figure 5-4. Definition of Sequence (uy (n*)) for N=4 
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Figure 5-6. 


Response of DFT FILTER to 
a Periodic Input for N=4 
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Figure 5-9. Availability Time of Convolution of Sequence 
Си (т)) for N=8 
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АРВЕМРОТ Хажд. 


CORRELATION IN LIMIT CYCLES 


E Introduction 

Limit cycles occur in digital filters as a consequence 
of the non-linearity introduced by the requirement to quan- 
Eze the results of multiplication. While analyzing limit 
cycles, Parker and Hess [52] showed that the quantization 
error sources exhibit limit cycles which are sinusoidal in 
nature but have a discontinuity if the principal magnitude 
is greater than the maximum quantization step size. For 
example, Figure A-la shows a limit cycle where the principal 
magnitude of the sinusoid is A while the maximum quantization 
error is M. Figure A-1b shows a iimit cycle where the princi- 
pal magnitude of the sinusoid is B but B is less than M so 
ao discontinuity exists. 

Returning to Figure A-la, notice that the quantization 
error departs from a pure sinusoid at most once,always re- 
turning at the next discontinuity. Now if A had been larger 
than 3M, the next discontinuity would have to move the quan- 
tization error farther Trom the pure sinusoid as Shown іп 
Figure A-1c. Іп this case the next two discontinuities make 
the quantization error coincide with the original sinusoid 
again. 

From the previous discussion, define a тїї ey ele Veen 


discontinuities which cause the quantization error to be at 
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most once removed from the sinusoid (Figure A-la) as a case 
l limit cycle. If the quantization error is twice removed, 
Bis a case 2 limit cycle. A case k limit cycle is then К 
Ж ИШЕЕ гешпоуей from a pure sinusoid; #f no discontinuities 
exist, this is case 0. 

Now consider two quantization error limit cycles in the 
same filter. It was shown by Parker and Hess that these 
will have the same fundamental frequency but one may be 
shifted from the other by a phase 6, as x, (n) and X» (n) ia 
Figures (A-la and b). When considering two error sources at 
a time, define the limit cycle condition as case К-2 where 
one error source Ís in a case k limit cycle and the other is 
in a case £ limit cycle. 

It is the object of this apnendix to show that aquantiza- 
tion error sources are not uncorrelated, especially in the 
limit cycle condition. The correlation coefficient for limit 
cycles of the 1-0 case, and the 2-0 and 0-0 cases are con- 
sidered, assuming analog rather than digital functions for 


ease of notation. 


ША Correlation in a 1-0 Limit Cycle 
Let the limit cycle functions for a 1-0 case limit cycle 


be given by the analog equivalent of Figures A-la and 4-1b 


and let 
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A sin wt CeO ey о ee {t,} 
or t € {(1-$)/w,(m+$)/w)} = {t} 


3 
xi (t) = ог © є {(2т-ф)/шв,2т/ш} = lto? 
Sin wt = eM te (9/0, =%) - {to} 
A sin ot + 2M Бе 1(п+ф)/о,(2п-ф)/о) = {64} 
x,(t) = B sin (wt + 9) 
where 


mit} represents the associated interval 

{-M,M} is the region of permissible values of x, (t) 
М <А < ЗМ 

0«B«M 

$ 


and 0 is an arbitrary phase shift. 


sin l(M/A), (i.e. sin l(1/3) < 6 « 1/2) 


The correlation between these two functions is defined by 


со 


E[x, (t£) x, (t)] A STE 


р 
12 010, 


st 
a 


со оо 172 
| Ј хубогрбае ухае) p(t) at} 


2T 
| (9 
s x, (t)x, (t)s pdt 
2т/ш 5 a 2т/ш 2 К 1/22 
2 Xx, (t)5 dt Xx) (t) 5 dt] 
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where folt) 15 the distribution of a random variable 6, 
uniformly distributed over one period of the signals. 


Performing the integrations yields 


1 
(r,6) = Іш - 8r (1-77) ^3cos 6 
M 1 1 
2 —- m —> 
[T - 16тг(1-г2)2 + 8mr^cos Kor ФО 71° 
тв г = “M/A, 5 ERA T 


3. Correlation in 2-0 and 0-0 case Limit Cycles 
When the same procedure is carried out for the 2-0 case 


limit cycle, the resulting correlation coefficient is 


1 1 

5 г г = 2 ~ i + РА 2 > 
(r,0) = ik iG eae a сов 8 
РОУ i 
[n^ 167^ (2n- sin") r-3sin"53r)-1607( (1-77)? 


NS 
aor) SE 


Bor the 0-0 case, there is no dependence on the ratio г = M/A 
since then A « M and no discontinuities occur in x4 (t). In 


this case the correlation coefficient is 
01248) - cos 6 


The correlation coefficient normalized by cos 9 is plotted 
in Figure A-2 as a function of a = 1/r = A/M for the 0-0, 


ED, and 2-0 cases. 
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APPENDIX B 


EXPERIMENTS IN ERROR CROSS-CORRELATION 


1. For Unquantized Signal Variable 


Нит co gain an insight into the nature of the 


, 


correlation between quantization error sources in digital 
weers, a simulation of the process of multiplying a single 
па! variable by two different coefficients and rounding 
the products was performed. A model of the simulation is 
shown in Figure B-1, where (u) represents a full precision, 
approximately uncorrelated Gaussian pseudo-random input with 
Zero mean and unity variance. 

The object is to multiply the input signal variable by 
wo coefficients, Az and Ал, and to find the erronos еј апа 
e»; by rounding the products (indicated by the operation Qu 


The correlation coefficient is computed by the definition 


E ou 
p = U . 
12 02 oc 
ш та 
where E[ - ] indicates “expected value" and А is the 
E 
Standard deviation of е, given by 
- m (В.2) 
ES d (ЕГе: 1) 


al 


mor zero meen variables. 
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To understand what occurs in the process, consider the 
error functions given in Figure B-2, where it is assumed 


that the ratio of the coefficients is 


2 


К = — = 


1 


г а 


(В.ч) 


and the As are both positive and rounding is assumed as the 
auantization method. 
The error е, is given by 


e, (u) ZEN u [A ulg EE k, (u)h «В 5957 


where h is the quantization step size and k, (u) is an integer 


determined oy 


k.(u)h Ee CoD 
ERU n a D B.6 
i SA < Uu € ———— t (В.ба) 


which implies 


A.U A.u 
i 1 і Ф} B) 
ELTE 3*2 PR 
ог 
Ru 1 
+ = 
k, (u) Бате a. 


zu 





where [ : Ју = greatest integer less than (+). It turns out 
then that the two functions е. (и) апа е (1) must coincide 
integer multiples of p = eh/A, = 3h/A,, since е. (и) 18 
periodic with period В/А- апа estu) has period Һ/А-. These 


periods coincide when 

Jh/A, = Kh/A, (B.7a) 
(where J and K are integers) which implies that 

EE = 273 ВТВ) 


by assumption. So the product e e, is periodic with period 


a 


р. 
When e,(u) is plotted against е. (и), the result is a 
diagram like Figures B-3(a through d) for the ratios given. 


If it is assumed that the errors are uniformly а СЕ 


Е еп -h/2 and h/2 and that the joint probability density 
is uniformly distributed along the solid lines of Figures B-3, 
mts possible to calculate the correlation coefficient 
meeording to Equation (B.1). 

For example: 


Consider Figure B-3c where R is given by (B.4). The 


joint probability density is given by 





1 Рог a continuous Gaussian input vith large enough 
Variance, this assumption is well justified. 
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p(ei,e5) 7 p(e5) · р(е, |е,) (B.8) 


where 
E 
p(e,) m (B.8a) 
lisle e -25)1) + óle MM 
3 jg БС ey7 l3 (26575) ] 
i §(e,-[3(2e,)])} T < ES <> 
ple, e,) = aU Ce, - L5 (26 5-h) 3) + 6 (е -СЪ(2е 2) (B.8b) 


+ 6 (е. -Г3(2е +) 1)) -7 S 55 <i 


+{5(еу-[5(2е„)1) + в(е-СЪС2е в) 1) 


4 Ша h h 
! 5 (с СЪ(2с 321) 127 2 < е; < Sm 
and ó(x) is the Dirac delta function. 
Then 
h h 
2 2 
т = 
Еее, ) / / е,е-р(е,,е.) de, de, 
EM 
2 2 
h h 
2 аи: 
= f е-р(е-) T е,р(е,]е„) de, de, 
h an 
ша 2 
h n me 
Ры» 2 Ц e, 
== / = де > + J ues de, T f ==? de, 
п 3h h _h 
2 
К ng B.9a) 
= - i 
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апа 


h 
e 2 
= 2 BE 
o. d e,p(e,) de, = ix (B.9b) 
та 
so that 
Е б 
Pio 1> (В.10) 


Performing the calculation for the rest of Figures B-3 


yields for 
R= 2: ру =- i (В.11а) 
ВеЗ: = = (B.11b) 
R = 3. ров = = (B.11c) 
R = 2: E i (Bere) 


It is suggested by Equations (B.11) that the rule Гог 
Ending the correlation coefficient is as follows: 
(1) Write Az and A, as integer multiples of their 


Bmantization step size, i.e. 


в. = Ja | (B.123) 


A, = Ка (ВЕТО) 


B _ _ 





~ 


Note that q is not necessarily equal родини 


(2) Reduce J and K until they are relatively prime 


So that the ratio of coefficients becomes 


A [K] 
6 К Т 
А g DI. 5p) 


where E means the product of the prime factors of m 


which are not prime factors of n. 


(3a) If UK]; and 71% are both odd, then 


signum [R] 


E NEST 
12 K oi J pK 


(ВА Да) 


иеге signum [R] is +1 for R positive, and -1 for R negative 
(К = 0 is a degenerate case). 
Әр)” ТІ ІК) is odd and 1710 is even or vice 


versa, then 


signum [R] 


1 
0 MEME (pom 
12 2 ТК Јој (к 


The computer simulation verifies these results. There; 


А. апа A, were made integer multiples of 0,01 between 0.01 


and 1.0. Ten thousand random numbers were fed into the 
simulation and the products were rounded to two decimal 


places. A11 data for the same ratio vere 222. and 


Ву he limits of the experiment, there was more data to 
average for low values of [J] „ and [X] ., thereby yielding 
better results.  Averaging wat necessarj due to internal 
computer rounding which aiso occurs and interferes with 
the simulated rounding. 


CM 





plotted in Figures B-4 through B-8. Only ratios greater 
than unity are considered and the data are divided into 
groups where the ratios are integers (LI) x - 1), multiples 
Df 3 (3l = 2 unless K is even), etc. (Data which appears 
оп а previous flgure are plotted in parentheses.) 

In each figure, the values of the correlation coefficient 


are also shown multiplied by (33,27. When this is done 


the points fall on the two functions 1/R or -1/2R. Thus 


VA [0] g/ LK 7 


2 


(LJ Jy) 9 or (B.15) 
-1/2R = го 5 
and therefore 
VARENNE SN 
р = or (ESO) 
-1/2121, [K] 


verifying Equations (B.14) except for the factor si num [R]. 
But it is obvious from Figures B-2 and B-3 that lof RS 


negative, Ay and As are of opposite sign and the sign of the 


correlation coefficient would be the opposite of what it 


1 апа A, had the same sign. (The experiment was 


performed for Ay and A, both positive.) 


would be if A 
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Eye Or eQueantized Signal Variable 

In Section 1, the distributions of the errors were 
assumed to be uniform. This is because the input had a full 
precision (essentially continuous) Gaussian distribution. 
When the input is quantized before the multiplication then 
the distributions and joint distributions become discretized 
and the results for the correlation coefficient are signifi- 
cantly different. А simulation of this case has not been 
run but the set of points in the (e, 85) plane where the 
Mont error probability density function has its support 
has been plotted for a signal quantization level of h = 1/16, 
a multiplier coefficient quantization level of q = 1/8. 
These are shown in Figures B-9. The support is indicated 
by a numeral in the grapn. 

The products (EJ)h = EJ/16 and (EK)h = EK/16 are the 
errors corresponding to the multipliers А, = Ја = Ј/8 апа 
A, = Ka = K/8. If it is assumed that every point of support 
ВЕ саца1 probability of occurrence, then the correlation 
coefficient ċan be simply calculated. By summing the pro- 
ducts of EJ/16 times EK/16 for all values of EJ and EK which 
provide support and dividing by the number of points pro- 
viding support, Еее, is found. Ther variance ot eden 
error can be found by summing the squares of the values of 
ЕК/16 (or EJ/16) for which there is support in that row (or 
column) and dividing by the number of rows (or columns) for 
which there is support. Note that it is always (де о о 


if there are n (>0) points of support in one row (or column) 


an 





then there are n points of support in any other row (or 
column) that contains support (as long as EK - 8 апа EK = -8 
are considered the same row). 

No attempt has been made to date to calculate the corre- 
lation coefficient for this case or to predict a formula 
for its calculation based only on knowledge of the relatively 
prime values of the multiplier coefficients, but this is 


suggested for future studies in this area. 


See conclusions 

a. Figures B-3 show the only possible values of е. апа 
es Miren can support their joint probability density function 
Tor the ratio given, no matter what distribution the input 
has. 

b. Equations (B.14) yield an accurate formula for the 
correlation coefficient between two sources of quantization 
error arising from the same signal variable lf the distribu- 
Mon of that signal variable implies uniform distribution 
of the error sources. 

c. Equations (B.14) do not reveal any dependence of 
the correlation coefficient on the input quantization level 
or the output quantization level. It is hypothesized that 
these dependencies do exist. In addition, there may be a 
dependence of the quantization level of the coefficients and 
the common prime factors of J and K. 

d. The patterning of the distribution support suggests 


ПИТТ су то the work cf Mulcahy [55]. He examined the 


correlation between the input to a multiplier and the 
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consequent error produced by quantization. His work could 
be extended and combined with the results depicted in 
Figures (B-9) to yield insight into the nature of the 


correlation between two errors. 
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Figure B-2. Error Functions for R= == 2 
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Figure B-3c. 
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Figure B-3d. 
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Figure B-9.18b. 
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Error Distribution for R = J/K 
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Figure B-9.19. Error Distribution for R = J/K = 5 
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Figure B-9.20. Error Distribution for R = J/K 
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Figure B-9.21. Error Distribution for R = J/K = 7 
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APPENDIX C 
VERIFICATION OF POSITIVITY OF FACTORS 
IN THE OUTPUT ERROR VARIANCE 

ime Introduction 

In Chapter IV several expressions for the variance of 
tie error at the output of a second order digital filter 
caused by quantization noise are derived. Since a variance 
is always a non-negative number, these expressions should be 
examined to ensure they yield non-negative results. 

All the expressions involve factors containing the 
КЕ Елее Ol the transfer function of фћенае рат ет as 
given in Chapter III. There are conditions on the coefficients 
in the characteristic funetion of the filter imposed by ine 
requirement that the filter be stable. The characteristic 
function of a second-order digital filter is given by 


22 мас р. Hor stability it is кгедашгед слаб =) < 5 < 1 


and -(1+b) <a < 1+b. There are no restrictions on the 
values which may be taken on by the coefficients c and e 


(or с! апд е!) since they do not effect the stability of 


the filter. 


e Positivity of the Factor, A 


Део поен ea appears in every expression for the 


Eutout error variance is given by 


2 


A * (1-b)[(14b)^ - a^] (C.1) 
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БИЕ | < о" < | it is clear that the factors (1-b) and (1+5) 
are both positive. The requirement that -(1*b) < а < (1%5) 
implies that la[ « |1+b| so that the factor in square brackets 
in (c.1) is positive. Therefore A is always positive inside 
the stable region. 

In addition to the positivity of A within the stable 
region of the (a,b) plane, it is obvious thet ^ is zero on 
the boundary of that region, since the boundary is defined 
by the equations b = 1, a = 1+b and a = -(1+b). Figure C-1 
shows the regions of positivity of ^. The points (a,b) = (0,-1) 
апа (0,1/3) are critical points of A. The point (a,b) = (0,-1) 
is a saddle point and (0,1/3) is a local maximum, where 


A = 32/27. At the origin, ô= 1. 


БОРО Ciney of the Factor, Сер) А 

Figure C-2 shows a plot! of the factor (1+b)/4. It is 
clear from the previous section that this factor is alweys 
positive inside the stable region since A > O and (1+b) > 0. 
Furthermore (1+b)/A has a critical point at (a,b) = (0,0) 


which turns out to be a local minimum. The value of (1*b)/^ 


at the origin is seen to be unity. 


КИМИ РОБ ТУ of the factor, p = eee Git )-2ace iA 


There are several ways to show that the factor, 


р = [(c^*e?^) (1*b)-2ace ]/A > 0. First of а11, note that it 





1 СВ Г. Souchon Federal German Navy, 
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is only necessary to show that the numerator is non-negative 
since, in Section 2 it was shown that А > 0 inside the stable 
region. 


a. One way to show that p » 0 15 to observe that 
(c^*e^) (14b) > 0 (С.2) 


so that it is only necessary to show that the largest value 
of 2ace is less than (“е^ iam». 

lf the product, ce, is positive then the largest 
дасе сап Бе 15 2се(1+Ь) з1псе а < 1+5. Letting с = ye with 


veo O forces ce = үе“ to be positive. Then 


(с2+е2) (1+) - 2ce(14b) > 0 (още, 

if and only if 
2, 

СЕЮ > оү (С.Д) 

anda (C.4) is always true. 
For ce negative, 2ace is maximum when a = -(1*b). 

Then to show that 

(с2+е2) (1+0) + 2се(1+Ы) > 0 (с.5) 


let c = -ye with y > 0. Then (C.5) is true whenever 


Ser 


а чн 





(1ғү2) Ec 


БЕ before. 

It can be noted here that the numerator of p is zero 
if c = e and a = 1+b or if c = -e and a = -(1+)). But then 
A is zero also and p ís a degenerate form. The numerator 
IESO zero far any a and 0 imside the s£2hl1e negian if 
све 0, 

b. Another way to show that р > 0 is to make use of 


section 3 by writing 
_ DAE 
р = (с *e^)[(14b)/4A] - 2асе/А. (с.6) 


For се positive 


2 2 


р > (сте -2ce)(1+b)/A > c te^-2ce = (zen) > 0 (Ce 


since a < (14b) and (1*b)/A Lp 


For ce negative 
22 2 
РВ ЕЕ ВС ТАБА > (cte) > 0 (С.8) 
since -а < (145). 
с. АНТА way to realize that p > 0 is to observe that 


КОО ШОК оао ос 1 e or e with the coefficient of the qua- 


Шабле бегм being (1+b)/A4 > 1 > 0. Fixing e and finding the 


328 





Born. where the slope of p with respect to c is zero, yields 
the value of c which make p minimum for the fixed e: 
. ae 
с* = „25 (089) 
Using с“ Бо define the locus of points where the 
mercial derivative of pP With respect to c is zero, substitute 
c* into p to form p* = p ы” TRAS 15 52111 а ашейт "с "Њу 
б 

e opening upward. Taking the derivative of p* with respect 


to e then yields the value of the parameter e which yields 


the minimum of the minima 


е# = 0 (C210) 


Which makes 


ее о (CRET) 
so that p is a universal minimum at c - e - 02 Теп р = 0 
unless A = 0. This case will be of no concern since when 
с = е = 0 the filter is degenerate and consists of, at most, 


a single delay element. It suffices to say that р > 0. 


BR Nature of the factor, pre 


The factor pte is given by 


Б (се?) (14) - гасе + сл 
А 


pte (ела 


Bo 





This factor is also a quadratic in c or e opening upward so 


the same techniques as used in Section 4c can be utilized. 


Taking the partial of pte with respect to e yields 


ac 


е* 


for the value of e which 
Enbstituting e* into ptc 


taking the derivative of 


-A(1+b) 


с“ = 
2[(14b)* - 


Using this value 


value of e 


а - 1) 
2 


(63) 
makes prc minimum for a fixed c. 
to form (ptc)* = (pte) | ox and 
(ptc)* with respect to c yields 


= ———- (С.11) 


of ec in (C.13) makes the critical 


CCS) 


Substituting (C.14) and (C.15) into pte yields 


Gare) -> = 


2 
= UN NP СОО ) 


4 


The minimum value of (C.16) in the stable region is 


at b = 0. Then 


pte E 
с=-1/2 
е=-а/2 
b=0 


-1/) COSTO) 





Although the factor ptc is not necessarily positive 
Bbself, the faetors in which it is used in Chapter IV are 


2(ptc) + 3 and 3pt2ct2 and by this analysis and Section 4: 
alpre + > 2.5 (CAME) 
and 
3р+2с+2 = 2(р+с)+р+2 > 1.5+р > 1.5 (6 19) 


псе р > 0. 


6. Other Factors of Interest 
The preceding sections have been concerned with factors 
which appear in the expressions for the output error variance 
for the uncorrelated case as listed in Table 4-1. There are 
factors in the general expressions for с 


Ау 
(4.48) for which knowledge of positivity may be useful in 


given by Equations 


Шеше studies. These w111 be discussed here. 
a. Іп Equations (4.48a and c) the factor multiplying 
f11 ШЕ аслошев со р іп Section 4+ of this appendix. This 


has been shown to be non-negative. 


b. The factor muitiplying f и Е ашавлоп (1. ЧВа) is 


да 


ӘӘ 222202220247 (156)-22(1-5)11/4 (С.20) 





This too is a quadratic in с or e opening upward. Taking 


phe partial ог г with respect to c yields 


с ае 
ck = Ed (C 2) 


ШЕ tune value of c which minimizes r for a fixed e. Defining 


rž = г | ах and differentiating r* with respect to e yields 
e* = 0 


so that c = e = 0 yields the minimum value of г. The factor 
ris unbounded on the boundary of the stable region unless 
р = -1 ог с = е = 0. Then it is an indeterminate form. 


ши бре factor ltatb multiplying 11 апа Er in Equation 


(4.484) is definitely non-negative by tne requirement that 
=(1+b) < a < 1+b in a stable filter. 


gee ine tactor multiplying f in (4.48c) can also be 


Е ща 
слот to бе non-negative, its minimum occuring for c = e = 0. 


Кле factor that multiplies f 2 іп (4,48с) has not been 


2 


Шшессеа сопр1Тесеіу, but it can be shown that it is a quadratic 
in c or e and that the coefficient multiplying either auadratic 
term is non-negative. It is expected that this term is 
non-negative. 


ШӘ те асбо (a) multiplying f Ма силе ом (41.48b) 


12 


can of course be positive or negative. 


EIN. 3ctcom multiplying f.. in (4,4588) can also be 


12 
positive or negative, since the coefficient which multiplies 


the quadratic term c or e can be positive or negative. 
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Fig. C-l. Sten of A = (1-b)[(1#b)*-a°] in the (a,b)-plane 
re ee е 
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APPENDIX D 


EXPERIMENTS IN ERROR PROPAGATION 
IN SELECTED DIGITAL FILTERS 


IEA TOUS cion 

Chapter IV presents theoretical formulae for the predic- 
tion of the variance of the output error caused by Quantize- 
eo or products in a finite precision digital filter. BMese 
formulae involve an infinite summation which must be checked 
for convergence. There is no doubt that the summation con- 
verges since it is assumed that the filter is stable and the 
entries of the correlation matrices are bounded. The Question 
is whether in the limit this summation yields the variance 
пе си вас еггог. But apriori sbatisties of the errons 
are not known, particularly with respect to correlation. 

In order to gain some insight into the processes involved, 
E filter were performed on an 
IBM 360. The objective was to test the convergence of 


computer simulations of a TM 


НАИШАО 21) as it applies to Equation (4.48b) for various 


inputs and sets of filter coefficients. 


2. Experimental Design 
Се Filter Simulation 
FOr E simulation, the canonice realization TMy¢ 
with d' = 0 was chosen because it does not involve quantiza- 
БИРИ ји the output equation, it has the simplest 


expres.ion for the output error variance, and it does not 
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пете any rearrangement of coefficients from those given 
in the transfer function. The state equaticns and the trans- 


fer function for this fliter are given by: 


x, (n) - -ax, (n-1) + c'u(n-1) (D.la) 
хо (п) = -bx,(n-1) * e'u(n-1) (DOES 
ие D (D.1c) 


and 


(ct + Mc 


l +t ame + bz 


H(z) (D.2) 


^ 
( 


A subroutine which performed two algorithms for the 
filter was designed. One algorithm utilizes double precision 
arithmetic operations and corresponds to Figure 3-42, the 
other corresponds to Figure 4-la with the errors being 
generated by a quantization subroutine which rounds numbers 
to one decimal place. 

роса с - 5 2305 Collected 

The difference between the signal variables of the 

two filter algorithms yields the accumulated errors Av(n) 


and y(n). The required statistics of these errors are 


Ше 
г = 


2 кш 
с (Ау) “с 


П 52) 


Гау(1)12- 15 Е дуа) ] (D.3a) 
J izl 


А 





1 


yGOy' G) - [y y (1) 1 


T 
(D.3b) 


И са ча 
и са 


у(1) (1 
T 


МУ) = 


<< 
= | 


Pel 


^ 
и = 
pa 


where M is the number of iterations performed. (Writing the 
statistics as being functionally dependent on Av and y indi- 
cates the source of the statistic.) The value of M for all 
ruins МОБ "9900. 

The quantization errors were preserved and all auto- 


and cross-correlation functions (as defined in Appendix E) 


for these errors were calculated up to Ik] 40 i.e. 


и са а 


ЕР 
Rig 6k) = e (kde (1) -40 < k < 40 (D.4) 


ізі 
where a, В 7» a,b,c' and/or e!. 

Furthermore, the quantization errors were summed as 
appropriate to form the composite errors e4 (n) and e,(n) as 


shown in Figure (4-1b) and given by 
e,(n) = e,(n) t e,,(n) (Пера) 
Eu nct e m) GN 
The auto- and cross-correlations of these errors 
were also calculated as in (D.4) except a,8 = 1 and/or 2. 


With these statistics available it is possible to 


form 
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^ m um m US t 

Hee) = Rago) + E A R (1) * ZO[ARQ)]' m=0,1,...40 
1=1 i=] 

| (D499 


а partial sum of the right-hand side 


where F (e) is the m 
of Equation (4.21) and depends onthe errors e. (Note that 
Е (е) contains only R,(0).) From (D.3b) it is also possible 


Lo form 
- E 
Е(у) = Л, (у) ~ АЛ (у)А (0.8) 


which corresponds to the left-hand side of (4.21). Then 
the convergence of F (e) to F(y) can be checked. 

For the output error variance the quantities с, (у) 
апа с, (е) меге calculated byosubstituting (D.8) апа (р.7) 
respectively into Equation (4.48b). The convergence of 
сд (е) Бо o£ (Av) eanscnensbesseen. 
| The auto- and cross-correlations of the input and 
the states of the low precision algorithm and of the state 
errors were also calculated and plotted. 

с. Filter Coefficients 

Three sets of filter coefficients were chosen to 
check the effect of changes in pole/zero locations on the 
error processes. These three sets are given in Table D-la, 
and the corresponding filters are labelled 0, 1 and 2 for 
bookkeeping purposes. Also included in Table D-1 are the 


resulting parameters of the infinite precision filter for 


the given coefficients where 


ES 


г is the damping coefficient 
w is the undamped natural frequency 


YT is the exponent of the envelope cver the 
unit pulse response 


ш is the frequency of the unit pulse response 
К, is a phasing constant given by Equations (».12) 
d. inputs 
Three different input conditions were used to derive 

the filters, a zero mean, unity standard deviation, essen- 
tially uncorrelated Gaussian pseudo-random one, and two step 
fünetionsşy one vith magnitudewl.0 and thesethems0.>b. Тһе 
random input was chosen to check the shape of the state and 
state error correlations. The step inputs were chosen to 
check the response to a highly correlated input. Since 5115 
can also be viewed — deterministic input, the output error 


variance should be interpreted as the mean squared errcr. 


-3. Results for Random Input 
The results of the simulation for the random input are 
shown in Table D-2 where the percent normalized differences 
are shown by 


07% - 100[o^ (e) - о, (Ау)/а, (ду) (D.9) 


2 
Ау 


the predicted value for the output error variance given by 


Scu) headed by "gp by Gre oom’. indicates 


Equation (4.66f) where it is assumed that the error sources 





аге uncorrelated and uniformly distributed so that 


д? = Һ2/12 - (0.1)2/12 = 0.08333... (0.10) 


The percent normalized difference for this theoretical value 
ENS given "m$ Column (12). 

All these data tend to confirm the convergence to 
the correct value. The most apparent result is that the 
rate of convergence depends very heaviiy on the vaiue of 
the damping coefficient, с. This is shown in Column (10) 
е) 
converged to and stayed within 0.1% of the actual value. 


where the value of mis given for the term when o 


Thus for filter 2 the lower damping coefficient required 
that more terms be summed to get within 0.12. It should be 
stated here that there appeared to be some divergence be- 
tween term 20 and term 40 for filter 2. This divergence 
was slow and seemed to be peaking somewhere near term 10. 
Another fact to be inferred@from Tables D-2 is whet 
the very slight shift in zero locations between filters 0 
and 1 resulted in a detectable change in the output error 
variance. This change is more obvious for the step input 
to be discussed in the next section. This change is not 
predicted by Column (11) since the "uncorrelated" assumption 
yields formulae which do not depend оп с! апд е!. 
miowter indication. that оне евтогэжаге пох unc omre— 
lated is demonstrated by Figures D-2 through D-22 which show 


the correlation functions of the states and state errors. 


ES 





(Recalling that for this filter the output is equal to state 
X4 ; figures which include tnis state equivalently include 

the output and figures which include the state error Уу 
equivalently include Av.) It can be seen from the autocorre- 
lation of x, in Figures D-2a, D-3a, and D-4a that the conver= 
gence of the output autocorrelation follows the same pattern 
as the convergence of the formula for the output error vari- 
ance. That is, in Figures D-2a and D-3a the magnitude of the 


autocorrelation of the output enters and remains in the 


region +0.075 at k = +3, whereas in Figure D-4a this condi- 
респ осопгв ас К - 213. The region +0.075 is also the region 
in which non-zero values of the autocorrelation of the input 
(shown in Figure D-1) occur for k # 0. This is a measure 

oz new close to uncorrelated the input is. 

These Figures also show the interesting fact that 
the accumulated state error correlations exhibit the same 
pattern as the finite precision state correlations which, 
by reference to Appendix F, can be seen to exhibit the same 
pattern as the theoretical infinite precision correlations. 
This is surprising since the quantization error sources and 
the composite errors еј апа е, did not exhibit any discern- 
ible pattern in their correlations (and thus were not inclu- 
ded here). There must then be some pattern in the error 
correlations which would result in the patterns of the state 
error correlations when properly combined. This combination 


might be found by examining the equation 
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Ry (Kk) B E(y(n-k)y" (n) ] 


n-k n 
а - А" Ан ебме (ее ул Ј 
1<1 1-1 


n-k n Ж 
= унунун 
151 j=l 


Me 91. (D.11) 
There are two more points to be made with respect 
Ene resulcs for the random input case, both involvime the 
Stevlsoles Of the quantization error sources. Fbst, ít 
was observed that the variance of all the error sources 
was very close to һ2/12 as given in Equation (D.10). This 
suggests that the errors probably were distributed in accor- 
dance with the assumption of a uniformly distribution between 
-h/2 and h/2. The mean of all the sources was zero which 
also agrees with this assumption. 
second, when the correlation coefficient for the 
errors ес! (п) апа e, 1 (1) was examined” it was found that it 
did not agree with that suggested by Appendix B. This con- 
firms the hypothesis of that appendix since the formulae 


given there did not hold when the input was quantized at 


about the same precision as tne filter variables. 


oe Пе о те о tioning of the multipliers c' and 
e! in Figure 4-1 is the same as that of the two multipliers 
in the experiment of Appendix В. 





4. Results for Step Input 

When the input to the filter was a step function of 
magnitude 1.0 or 0.5 some significantly different results 
were achieved. Table D-3 lists these results for mean 
squared errors rather than variances. 

Column (11), headed by 'hv^(e) by (4.91b)", indicates 
the predicted value for the mean squared output error given 
by Equation (4.91b) if the correlations are assumed to be 
cemevanc and egual tothe a posteriori vare ас к= о. 

The percent normalized differences given are calculated 
in the same way as in (D.9). 

Again the data confirm the convergence to the correct 
value. The rate of convergence also depends here on the 
damping coefficient, 5, but compared to Table D-2, more 
terms are necessary because of the highly correlated nature 
of the errors. This is shown in Column (10) by the value 
of m required for convergence to within 0.1%. 

Шорасеношас the small shaft in zero locations between 
filters O and 1 caused a significant difference in the 
resultant mean squared output error. 

The plots of the correlations for the step inputs were 
nearly constant for the states, state errors, quantization 
error sources and composite error sources. There was a 
tendency for these correlations to be less than the center 
value at |k| = 40, but this deviation from an exactly constant 


correlation is to be expected for e real situation since no 
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real power spectrum can be the ideal Dirac delta function. 
The deviation from a constant value was less than 1.0% at 
[k| = 4o. 

| thesamean squaredeenrproreof theserror-seurces was net, of 
course, equal to һ2/12 for this case since this is more а 
deterministic case than a probabilistic one. Rather, Lesi s 
most likely that the states were very nearly constant after 
the transients died out and that the error sources behaved 
similarly. Thus the only limit cycle present was a "sticking" 
effect, that is, a constant difference between the "infinite 
precision" model and the low precision model. Тһе mean 
squared error of the error sources would then depend upon 


the value at which the states "stick". 


5. Suggested Changes in the Experimental Design“ 
аА сег Structure 
It is recommended that filter algorithms be developed 
for the rest of the canonic realizations of Chapter III in 
order to compare the results for different realizations. 
b. Quantization Subroutine 
Some error was introduced into tne experiment by 
having a quantization subroutine which attempted to round 
to some number of decimal points while the program was run 


on a machine that operates in binary. It is a simple matter 


“The program used for this experiment is available from 
E" c ar ker, Electrical Engineering Department, Naval 
os beEraduate School, Monterey, California. 


p 





to rewrite the subroutine to round to some number of binary 
points without effecting the basic principle of the experi- 
ment. The error experienced with the present subroutine is 
probably in the fifth or sixth decimal place which can be 
significant for a case like that of filter 2 with a unity 
step input where the mean Squared error had its first signi- 
ficant digit in the fifth decimal place. Additionally, sub- 
routines which perform truncation and sign-magnitude trunca- 
tion шема be useful. 
с. Statistics Gathering 

The calculation of variances in this experiment was 
valid in that no variances were computed until the transients 
of the filter nad decayed (after forty iterations). But 
the correlations caiculated did make use of errors generatec 
during the transient period. Although 3000 more iterations 
were performed, thus making the effects of the transient 
errors small, an even more accurate experiment would result 


if these errors were not included. 
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APPENDIX E 


SPECDRAL THEORY ТМ DIGITAL FILTERS 


Porro dee tion 

In fields employing analog systems, a great deal of use 
has been made of spectral theory in dealing with stochastic 
processes. In digital signal preeesséng,- although.all the 
necessary tools for using spectral theory are well-defined 
[5,56], it has not been relied upon very heavily in this field. 
In this appendix, these tools are presented and will be 
utilized in Appendix F to gain knowledge about the spectral 


properties of the signal quantities in a digital filter. 


2. Synopsis of Spectral Theory іп the Z-Domain 
а. Power Spectrum of a Single Random Process 
Given a wide sense stationary random process, x(n), 
its autocorrelation function is defined as: 


4 


R_(k) Á Elx(n-k)x(n)] k,n integers (EL) 


Since this is a'function defined on a discrete domain 


it is possible to define its two-sided Z-transform as 


A NE: -k 
S (z) 8 Z[R,GOJ = ү UO в 2) 


2442 is called the power spectrum of the random process. 


When normalized S (z) is called the power spectral density. 
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The reason for this terminology is seen when one examines 
the relationship between R, (0) and the integral of S (z) 
around the unit circle in the z-plane. 

Since В (К) is the inverse Z-transform of S (z) it 


рат be written as 
Be): к-1 
ROO = aig j s, COz az (E.3) 
(All closed contour integration in this appendix is counter- 
clockwise on the unit circle.) 


Then the power in x(n) is given by 


E[x“(n)] 


— dz 
R, (0) = 3114 5,00 


ТІ 
ДІ; JOT 
A E 5 (е “)аш (5.4) 


which is the mean square value of the random process. 

| Other properties of the autocorrelation functions 
are that it is an even function and that the value at the 
origin is greater than or equal to the magnitude of the 


function for any other value of the argument, i.e. 
R,(0) > |R,(k) | (Е.5) 


Properties of the power spectrum include: 
(1) It is real and non-negative on the unit circle 
for real random processes. This follows since К (К) is real 


and even. 
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(C —Xt Ts ameven fumetlon of” w on the unit oebtrole 
Cz Е popup" In fact, it has even symmetry with respect to 


inversion of the argument, i.e. 
S (z) = S (1/2) (E.6) 


This also follows from the evenness of R (x). 
b. Cross Power Spectrum of Two Random Processes 
Given two random processes x(n) and y(n) which are 
mutually wide sense stationary, their cross-correlation 
mMaction is glven by 


Ryy(k) & ELx(n-k)y(n) J (E.7) 


The power spectrum of this function is also given 
Dy its Z-transform 


со 


-К 
= № 
Зху(2) A NO (E.8) 
The cross correlation function is not necessarily 
even, but it obeys the property that taking the cross corre- 
lation of the random processes in the opposite order has the 


Siibeet Of flippine the Minetlon@abeut the vertrealeaxis, i.e. 


Rx, GO - ER (E.9) 
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The magnitude properties of the cross correlation 


еке; 
(ву ПО | < 588, (0) + R, (0) J (Е.10) 
апа 
2 
[ROO | < R,(0)R, (0) (E.11) 


Property (E.9) suggests that the cross power spectrum 
is not necessarily real or even on the unit circle, but that 


same property implies that 


S ytz) = 5 уу(1/2) (E.12) 


but 


en 2 Sky 1/2) 
necessarily. 
c. Power Transfer in a Digital Filter 
(1) For One Input and One Output 
Let u(n) and x(n) be two random processes which 
are related by 


со 


x(n) = E u(n-m)h(m) (Е.13) 


m= EXC 
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ог 


X02) = Не =) (Е.14) 


That is, u(n) is the input and x(n) the output 
of a digital filter with unit pulse response, h(n), and 
transfer function H(z). Then the autocorrelation of x(n) 


is given by 


К (К) Е|х(п-к)х(п) | 


Ef E u(n-k-m)h(m) E u(n-a)h(q)) 


m= —oo а==% 
= у г R, (k+m-q )h(m)h(q) (E.15) 
== 128) qu 
and the power spectrum is 
и) = Е R eus 
X DER X 
=> i Y R (k-g+m)h(m)h(a)z”“ 
== == СО m= =% q--9 
со а со = È 
= у US Y h(m)z у п(а)29 
= S, (2)H(2)H(1/2) (ЕСЕ! 


The product H(z)H(1/z) is called the power transfer function 


for obvious reasons. 
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(2) Between Input and Output 
The cross-correlation between the input and 


output of the previous section is given by 


E[x(n-k)u(n) ] 


FM) 


| 
62 
ren 
са 


u(n-k-m)h(m)u(n)] 


1 


L R,(ktm)h(m) ET) 


m= =co 


and the cross-power spectrum is 


8 


Sure) 


-k 
= R (ktm) h(m) z 


1 
и са 8 


~ 
| 
8 

23 


ИТ ~ 


= СО 


со 


К (ра? у h(m)z" 


--СО = со 


| 
1M 6 


р 


S, (2)H(1/2) = S, (1/2)H(1/2) (E.18) 


When taken іп the opposite order 


co 


В х СК) у; В (к-т)а(т) (Е 19) 


= 


| 


S.) 


m 5 (2) Н(г) - 5 (1/2) (Е. 80) 


(3) For One Input and Two Outputs 
Let u(n) be a random process used as the input 


ОО ШОШО а filter with two outputs (or equivalently, to 
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two digital filters) and let the two outputs, x(n) and у(п) 
be related to u(n) by the transfer functions, G(z) end H(z), 
respectively. Then the cross-correlation between x(n) and 


y(n) 


ка 


m E[x(n-k)y(n)] 


E{ 2 u(n-k-m)g(m) E u(n-q)h(a)} 
m= =% q ==% 


= 2 R,(k+m-q)8(m)n(q) 


and the cross power spectrum is 


со со со 


5. 402) E L L R (k*m-a)g(m)h(q)z ^ 
xy kz-o m= =% а=-9% u 


L В (р)? у (т) = 
p ==% m--99 q 


| 528 


h(q)z 3 


S, (z)G(1/z)H(z) 


When taken in the opposite order 


со со 


А X OO di. os В (к-тта) в() к(а) 


uz) 


5 (2) С(г)Н(1/г) = S ,(172)G(z)H(1/z) 


| 


Seel) Az) 


xy 
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APPENDIX F 


SIGNAL CORRELATIONS FOR A TH C FILTER 
E introduction 


t 
00 


form oT a second order digital filter as described in Chapter 


The computer simulation in Appendix D utilized the TM 


TTI. ~The program-produeced plobs-of the-autocorrelation and 
cross=correlation functions of various signal and error 
quantities involved in the processing of uncorrelated and 
correlated inputs. Іп this appendix, the spectral theory of 
Appendix E will be used to derive the expected results of 
the plots had they been performed for an infinite precision 
filter. Deviations from these predictions are discussed in 


Appendix D. 


: ; Coe 
2. Unit Pulse Responses in a ТМ о Filter 


БЕ сла со Chapter III, the state equations for a TH C 


filter are 


x4 (n) = -ax, (n-1) + x, (n-1) + с!и(п-1) вета) 
X5(n) - -bx, (n-1) * e'u(n-1) (Кр) 
пл - x4 (n) + а'и(и-1) (Е.1с) 


Letting H4 (z) and HA(z) be the transfer functions from 


the input to x4 (n) and xin) respectively then 


ENT 





1 


E а т 1271 


H, (2) Б = (Е.2а) 
l + az 7 Be 
1 no t ei 
H (z) = C RES ae (F.2b) 


1 + az + bz 


Note that for 4! 03 H, (2) = H(z), the transfer-function 
of the fidber, but femad' - +, H, (2) is as given in (F,22) 
рим, the transfer function to the output is 


AA Sen кесі 


H(z) =z Е - 27 % Н1(2)  (F.3) 
1 Жай» + biz 


Now consider the Z-transforms of the two functions 


ә ав BnT n > 0 
g (n) = (F. 4a) 
0 n < 0 
ВОИ ТИ gnT п> 0 
в (п) = (F.4b) 
0 п < 0 
which аге 
E E 
M Ка 
С. (2) = (F,5a) 
l 1 - 2e eos gTz = 


e Ylsin grz 1 


а, (2) = a =e Se Е 56) 
- Т2 De ROS BT z te eyT,-e 
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Any transfer function of the form 








(c Ep TRE 
G(z) = mw - 2 (Раб) 
Ра? 
can be written as 
-1 
= + 
С (ра), c [G (z) ко, (2) ]2 (57) 
by writing 
с 
2 =l E 
(с са) а 1 e, (1 + o )2 = 
G(z) = CN -> = HA =” = с) [6., (2)+К6 (2) ]z n 
d puc 507 1 + аг Tob = 
с (1 - - "ub c. gTz "i к ке Тізіп a ғ. 
- Ент ORs eg MU дано ла ӨЛ 
= YN . Y / 
E 2e" YT cos BTz Е meii. : 
Equating coefficients it is found that 
m l 
np um b (F.9a) 
ВТ = msan E ) (F.9b) 
с 
2 + е с ВТ 2c, - ac 
с. 2 1 
В = — (F90) 
e sin BT PEE 231/2 AVR 
ШІ d 


Here 8 is the discrete response frequency o, апа у = zw 


а 
where t is the damping coefficient and o is the undamped 


ват Trequency of the filter. 
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Taking the inverse Z-transform of G(z) then 


g(n) = c,lg, (n-1) + kg,(n-1)] 
o et T tees B(n-1)T + ksin 8(n-1)T] nc 
0 15-50 
GEO ) 


For the transfer functions of Equations (F.2), the unit 


pulse response of the states of the тм о fever are 


euasit) ces B(n-1)T + xsin Во ccc 
h, (n) = 
0 Шы 0 
(Қаша)? 
етет Ќ0-1)7[сов Вп-1)7 + коза Вп-1)171 n» 1 
h (n) 
0 n 0 
(ЕО) 
where 
де! - ас! 
к, = кошы аа. аы (Е.12а) 
2е! (ъ - а)1/2Уъ 
ы 
k = ае! - гђе! (F.12b) 


2 pet(b = к ln 
EIN 


and y and ß are given by (F.9a and b). 
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3. State Response to an Uncorrelated Input 
When the input, u(n), to a digital filter is a wide 
sense stationary uncorrelated random process its autocorrela- 


tcron function 2S "given by 


No k = 0 
в, (к) = Elu(n-k)u(n)] = (F.13) 
0 k # 0 
where No is the mean square power in u(n). Its power 
spectrum is then 
i k 
= г 1 = F = 
5 (2) ZUR (k) J ae No CR IA) 


Referring to Appendix E, the power spectra for X4 (n) 


t 


00 Aten are 


and x(n), the states of the TM 


т - S, (2)H, (2)H, (1/2) = NH, (23H, (1/2) "“,15) 
I - S Cz)H, (172)H, (z) - МОН. (1/2)Н. (2) (F.16) 
"ue - S, (2)H, (2) = Мон. (=) (Т) 


Бог 4а! + 0, the output power spectrum is given by (E: 15) 
with i = 1 since then v(n) - x4 (n). This is aiso the cross 
Spectrum between the output and x(n). The eross spectrum 
between the input and output is (F.17) with i = 1 and the 


cross spectrum between the output and хә(п) is given by (P. 16). 
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The auto- and cross-correlations involved are by defini- 
tion the inverse Z-transform of the spectra. But observe 
that each of the spectra is the produet of two funetions in 
the z-domain. By the convolution theorem this corresponds 
to the convolution of the inverse Z-transforms of the two 
functions being multiplied. If the inverse Z-transform of 
H, (z) is h, (k), then the inverse Z-transform of H, (1/z) is 


h, (-k) = орао 


В 
МН (К) 


E _ 
2° (Hy (z)H,(2/z)] = hy Ck) *h, (-k) 


8 


Е ав (сет) = 


(m-k) 
(EE) 


Шам 


Zum 


which is the correlation of в; (К) with h, (k). Since Equations 
(F.11) require that h,(k) be zero for k < 0, (F.18) can be 


written as 


eo 


X h,(m)h,(m-k) ЕО 
m-l * J 5 
ER, (k) = со 
o + zo h.(m)h,(m-k) k > 0 
m-k*l ^ Y 


8 


h. (m+1)h, (m-k+1) k < 0 
m=0 + Ј т 


8 1 с 
© 


(Е.19) 


ИМ 
O 


Ў h, (mF1)h, (n-k*1) k 
E J 
m=k 
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Considering the functions of Equations (F.11), then 


Noh, (m+1)h,(m-kt1) = Kyo eos BmT + k,sin BmT ] 


[cos &(m-k)T + k,sin B(m-k)T] 
(Е.20) 


Since the indices п+1 and m-ktl are at least unity for the 


values of m and k in (F.19), aná Ky is a constant involving 


N, and the scaling constants c' and e! in (F.11). 


0 
By a change of variable in the second summation of (F.19) 


the correlation function R, j (k) becomes 


(> 


со 
-=Y|k|T -2ymT : 
e Y | | ME [cos MT + k sin Se! Wc 


K 4 
1j m-0 


-[cos ß(m-k)T + k,sin B(m-k)T] 


DN m 


ос 
ет КИ 2 ee B{mtk)T * k,sin 8(m*k)T] 


K d 


ES m=0 


[Г Өс + k sin Вит ] Е 7 0 
qon) 


It can be shown by the use of trigonometric identities that 


E үрт 
К; (К) К је (6.005 ВКТ *,, in Boum CES» 
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where 


-2YT -2YT 
б (ПЕТА )(1-е cos2BT)+(k,+k,)e sin2BT к К 
E wr. m а ии. 
Dm EI gu Me 


Теге cos2BT+e 128 


CE) 


| а ) (1-e ^ совавт)- (1-кк.) ет? sin28m к, К. 
2 2yT GYT E 


1-2е cos2ßT+te ти 
(F.21) 


Мите thespluswsien for k-- 0 -andethe mimis- for k < 0. Notice 
that for i = j, Ra 4,00 is an even function of k and its maxi- 
mum occurs at k = 0, so that R, 4 (k) is a valid auto-correlation 
Sum воле 


The point to be made here is that К, (k) is bounded in 


magnitude by кр ^ ES Pēr i = j, ths poured 
is unambiguous. But for i # j the choice of the sign in the 
expression for Ui; must be the one which maximizes Viso None- 
theless, this Bound Пав the same form for every correlation 
function involving the states and output of the infinite 
precision тм, filter with а! = 0 when the input is an 
uncorrelated wide sense stationary random process. А similar 
analysis would show that this is true for any filter realiza- 
tion with the same characteristic function. The only excep- 
eionelssechat at9ethe origin, К, (к), the output autocorrelation 
will have an additional constant term when d (or d!) = 1. 


This occurs because а = 1l implies a transfer function of 


unity between the input and output in parallel with the 
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rest of the transfer function. Then tne input autocorrela- 
tion pulse at k 0 is transferred directly to the output. 


MO snow this, let 


ee = H(z) = 1 + G(z) (Е.25) 


waierewaz)_ ds in the formeof (F.6). Then 
H(z)H(1/z) = 1 + G(z) + G(1/z) + G(z)G(1/z) (Е.26) 
and 
= 
В (к) = 55; 9 Н(2)Н(1/г)4г/2 E27) 


for an uncorrelated input so that 


со 


1 + 2g(0) + I g^(m) k-0 
m=0 
== (К) = | g(-k) * EX g(m)g(m-k) k«0 (F.28) 
o" m=0 
e(k) + £ g(m)g(m-k) ie Sf 


m=k 


исе (п) = О forn < 0. So fork 7 0 RR) is bounded 


by a constant times сет Оц а оне отп Сеге iS ап 


additional [1 * g(0) JN. term. 





ц. State Response to a Highly Correlated Input 

In the previous section an uncorrelated input was used 
gasdecermine the correlation Tumetisns of the states of the 
filter. The input power spectrum in that case was found to 
be a constant. At the other extreme, suppose the auto-corre- 
lation of the input is a constant. This is an extreme case 
because as noted in Appendix E the magnitude of the auto- 
correlation function at the origin is greater than or equal 
to the magnitude at any other value of the argument. It 
will be shown here that this kind of input yields correlation 
functions of the states which are also constant. 

dune Power Spectrum for a Constant Autocorrelation 

The first problem involved in finding the рөтег 

spectrum of a constant autocorrelation function is that a 
constant discrete function is not absolutely summable so that 
the application of the two-sided Z-transform may not be 
valid. This difficulty can be removed by letting the 


constant autocorrelation be given by” 


R,(k) = Tan ES a > 0 (F.29) 
о->0 


It is now possible to apply the two-sided Z-transform 


to 


Imhis method was suggested by Pror. Craig Comstock, Car. 
USNR, of the Mathematics Department, Naval Postgraduate School. 
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Да. ча 
R (k) =e mad (F.30) 
since this function is absolutely summable. 
The power spectrum of ко (к) is 
- k| -k 
54 (2) = 2[8°(к)] = z e^?lklz- 
u ц £ 
k=.0 
- К-К а кк 
= È ere 2” са 2: ess 7 
k=0 k=0 
-Q 
10 е г 
= — + - (ua) 
1 - е a 1 - е A5 
a a 


as long as z is in the annulus e`” < |z| < e“. Then by 


collecting terms 


AR 


3%(2) = 
а CIS S ER 


СЕЕ 29 


Letting Sz) be the analytic extension of the power spectrum 


and applying the inverse Z-transform, then by residue theory 


al а k dz 1 (са Се?) 
в с = ба? 
u 21, 40292 2 2j 0 ce m E ща 
сле k > 0 | 
а = elk! (Е.33) 
po k < 0 | 
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SOM ches is a Consistent transform pair. Notice that in 


the limit the power spectrum БИ (г) becomes the Dirac delta 


function at z = 1, since 


S (z) £ Z[R (9) = Z[1ám БХК) ] 
u u ap Ч 
= E lim e^ Ikl,-k = lim 2 Жы -k 
k-0 а>0 · с>0 k=0 
d со =] 
= Jim S (z) - = 6( 27-1) 
a0 M 0 otherwise 


This yields the transform pair 


ROO 7 S, (z) 


1 + 6(z-1) 


1 k = 0 
> 1 
0 k ¥ 0 


b. State Correlations 


It is now possible to find the correlations among 


Cie 


(Е. 


(F. 


Bo» 


CA) 

u 
c 
~“ 


350) 


the states and output for a constant input autocorrelation. 


In any second order digital filter, there are two transfer 


funetions between the input and the states and one between 


input and output. Let them be denoted by H(z), H, (z) and 


H (2) respectively. 


From Appendix E, the power spectra for the states 


and output are given by 
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814 C2 = SU UH, GZ22H, (2) (ED) 


where i,j = 1,2,v and 5, (2) is the power spectrum of the 


роде. Kor the case of the highly correlated input, let 


= о ss Е а 
S142) - Lip S,,(2) = ui 5402) болс 


and consider 5,22). The corresponding correlation function 


is then 


а A A а К dz ж 
К; (К) = 517 б Hi (72)H, (2)8, Cz)z. та (2.38) 
Now for any second order filter the product A, (172)R, (2) 


can be written as 


H, (1/z)H, (z) = Ki, (2) (Р. 39) 
j E 
e 


1 )(z-z 


(z-z,)(z-z,)(z-z 


where 24 апа 2. 


a polynomial in z. Then the integrand in (F.38) becomes 


are the poles of the filter and К. (2) > 


I(z) н(а)н (а)50( а) 


к, Са)(е 0 = са и 





(z-z,)(z-z5)(z-z] )(z-z5 )(z-e ?)(z-e^) 
(F.40) 
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For k > 0, it 15 possible to write the partial ехрап- 


Seon of (F.40) as 


ЕВ к) Вы в. (к) | CQ(k)  C.(k) 
T(z) т (ee) == М бы | - E - -1 b - Б с 


230607 











~Q а 
BOS Z=e Zt 


(Е.41) 


where С. (к) and C Ck) contain the factor (e “-e%) in their 


denominators while the other residues do not. Letting 
с. (к) С. (к) /(е “-е) (F.42) 


then by residue theory 


Ry (kK) [A, (k)+A,(k) ](e eS) + Ci(k) k > 0 


(F.43) 


so that 


m= lim DES = л с! (к) - Jim I(z)(z-e 9) 
2 део EU a0 z-e 


p eH, (e@)H, (e7%) 
v Ј 


(F.44) 


V 
© 


ii 


Hy (1)H, (1) к > 
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It can be shown that Equation (F.44) also holds for k < 0 


except that the last limit is 


-аК а - 0 ак -QK 
E (е H, (e ЈН (е ) + (е EC )H, (@)H,(0)} 


= H,(DDH,.(2) (F.45) 


as long as the filter has no pole at z = 0, 

The final result” is that if the input to a digital 
fuer nas“arconstant autocorrelation function equal te No» 
then the auto- and cross-correlations of the states and output 


mor Sven by 
Ry į (k) = Nol, (1)H, (1) (F.46) 
(where i,j = 1,2,v), which is also a constant correlation 


function. ЈЕ can also be shown that the cross-correlation 


between the input and one of the states or the output is 
со 
В (К) = М.Н, (1) = № : h, (n) (F.47) 
Notice that the cross-correlation functions are even func- 
tions in this case. Notice also that the quantity H, (1) is 


the sum (discrete integral) of the unit pulse response, hi (n), 


of the associated state or output. 


DOT 





APPENDIX G 


TRANSFORMING (m)"" ORDER COMPLEX FILTERS 
INTO (2m)"À ORDER REAL FILTERS 


In Chapter V, the DFT Filter is an algorithm which con- 
tains first order complex coefficient digital filter sections. 
These filter sections are purported to be digital resonators 
at certain frequencies around the unit circle in the Z-plane. 
But a resonator is expected to be a second order real coeffi- 
cient algorithm. It is shown here that the equivalence is 
justified. 

First consider the first order digital filter using 
complex arithmetic with one input and one output (see Figure 


G-1): 


x(n) ax(n-1) + bu(n) га. Ша) 


y (n) вн) T dulnn) (С. р) 


шестте нае u, x and y are “complex numbers. Let 


x(n) = x4 (n) + jx5(n) (G.2a) 
y(n) = y} (n) t jy, (n) | (G.2b) 
u(n) = u, (n) + Ја, (п) (G.2c) 


Боје 





c + да, (G.2d) 
р = bi + Jb, er 
с=с, + је, (пег) 
Q cma, није, (G.2g) 


where subscript (1) represents the real part of the complex 
number and subscript (2) is the imaginary part. Then Equa- 


tions (G.1) may be written 


x4 (n) + jx,(n) a, x (n-1) - а-х-(п-1) + JLa, x; (n-1) 


+ арх. (п-1) J 


+ 


Би. (а) - b,u,(n) + j[b,u, (n) + bou, (n) J 
(G. 3a) 


у. (п) + jy, (n) = сх. (п) - c„x,(n) 55 jle,x5(n) В сх. (п) ] 


+ 


d, u, (n) - аи, (п) + jld u, (n) + d5u4 (n) ] 
(G.3t) 


Equating the real and imaginary parts of Equations (G.3) and 


using state variable form yields: 
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x4 (n) а, -а, x4 (n-1) p WC u, (n) 
= + 

a) a, а. | | x5 (n-1) b, bi u, (n) 

(G.4a) 

y, (n) су -с- x4 (n) а. -4, u, (n) 
= + 

y(n) с C4 X, (n) а, а. u, (n) 

(6.545) 


These can be recognized as the equations for a second order 
digital filter using scalar arithmetic with two inputs and 
two outputs (see Figure G-2). Тһе characteristic equation 
of this second order system is: 

АЕ 


TY 2a,z + а + as к 22 - eR,lalz + Е 
The pole locations are therefore z = a and 2 = a* where * 
indicates complex conjugate. Note that if (a) is a real 
number then a, = 0 and a = a* and this second order filter 
is actually two uncoupled first order filters each with a 
pole et Zz = а = ај. 
The general result is that an ао order digital filter 
using complex arithmetic with p complex inputs and q complex 
outputs can be transformed into a (от) order digital filter 


EP rrvehmetiec with 2p scalar inputs and 2q scalar 


outputs. This is shown in the following development. 
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For the па 


х(п) 


y(n) 


where 


u(n) is 
y(n) is 
А is an 


B is an 


order complex filter we have the equations: 


Ax(n-1) * Bu(n) 


Cx(n) + Du(n) 


is an m x 1 column vector 


ap x 1 column vector 
a q x 1 colum veever 
m x m matrix 


m x p matrix 


C is a q x m matrix 


РБ тара тарна 


х(п) 


Let 


Rx, (n) be the real part of х (п) 


x, (n) үү 
x5 (n) Е сащ 
x62) | “m1 


(G.5a) 
(GWS ) 
12 кл 
22 PARI 
Е | Е 
mm | 
(6.6) 


Ix, (n) be the imaginary part of x (n) 


На. + Бе 


the real part of arj 


> 





E be the imaginary part of a 


еее. 


ij 


Let RA - (Ra, ,) IA = (Ta, ,) Rx = (Rx,) etc. 
i.e. RA is the matrix whose elements are Ra, у. Then equating 
real and imaginary parts in (G.5) and forming partitioned 


vectors and matrices as follows: 


Rx(n) ВА | -iA | | Rx(n-1) RB , -IB | | Ru(n) 
----- = ee: ————— 3E ro — == 
| 
Ix(n) TA ı RA | Ix(ne2) IB ! RB ||IuQ) 
(G.7a) 
Ry(n) RC | -1с | | ха) RD | -ID || Ru(n) 
---- = E ---- | | ------ T ETE — 
| 
Ту (п) О од тост ID ı RD ||Iu(n) 
(6.70) 


Another form of these equations is found by replacing each 
complex element in each complex vector and matrix with its 


equivalent scalar vector or scalar matrix defined as 


Rx, (n) 
х,(п) > а х (а) (6. 8а) 
1 Ix, (n) E 
Ras 221 
id и : A (G. 8b) 
ја + На + 
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Then Equations (G.5) become 
X. (n) А EE A Бас 1) ВоВ | IB 
ED 411! c o | 2 а EE n U Cn) 
| ---1---4-- +--- | |----- 
Ша 2. spec c MM EN Bo | 
ася 221, 222 5 Ша 
| ---.---к- = 
X,(n) ae Aso . X, (n-1) 1 . | . . 
um mm = = о = е ЕЕ о С 
= . | 2 + ° 
— 8 А A 
X (n) A. А X (n-1)] IB | | U (n 
= -ml! | mm 5 =m! ¡mp | Ше ) 
(G.9a) 
vl ~ | | Wa 1 Г | | | 
КИ бре Sm EET 15, |1,0)! 
ШЕР; ou i © Cae ENT cw ENDE EDI. 
el от | 
сш ч, К 
! U, (n) 
= ' tL —— . Ho. 
| у 
ань am am m = а аи өн» | ань ань ешь 
dle nn 
| i || 
к (а) e D | Сал X (n) Doi (Pap jj 95 02? 
| | | 
(G.9b) 


Equations G.7a and G.9a each represent 2m real state equations 


with 2p real inputs. 


2q real output equations with ?p real inputs. 


ET 


Eguatdons G/D and G.9b each represent 





х(п) 6 о y(n) 


Фо qe 
o 


O 


Figure G-1. First Order Complex Coefficient Digital Filter 


u(n) 





Figure 6-2. Second Order Real Coefficient Digital Filter 
Equivalent to Figure 6-1. 
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APPENDIX H 
COMPLEX INTEGRATION INVOLVING A 
COMPLEX DIRAC DELTA FUNCTION 
From the results of Appendix F it is possible to propose 


the following theorem. 


THEOREM: If H(z) is an analytic function on the unit circle, 
then if H(z) has only simple poles inside the unit 


circle 


6 Н(2) 6(z-1) = = 2rjH(1l) (H.1) 


where the integral is taken counter-clockwise around the unit 


circle and 6(z-1) is a Dirac delta function defined by 


= а со даи 
8(z-1) Ê lim EB. == р = (н.2) 
а+0. (2-е ^" )(z-e ) 0 2 1 
PROCE: Consider the integral 
слави о 
Есе а. o (8.3) 


(Се uS) 


if ОЛЫ NM sumple poles Inside the unit circle, then by 


residue theory 


Т(а) = 2пј(е се“) R. * 2njH(e" ^) (Н.Д) 


1. 


Ub 


i=l 
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th 


where R, is the residue of H(z) at the 1 pole 


1 - а 
(zee )(2-е ) 
of H(z) inside the unit circle. Notice that R, does not 


i 
have the factor (а Чес) in its denominator unless H(z) 
has & pole at C To avoid any conflict assume that the 
pole with the largest modulus inside the unit circle has 


|z| = b, and restrict a such that 0 < а < -1п b. Under 


this condition 


lim I(a) = 21jH(1) (Н.5) 
Q0 


On the other hand, by (H.3) 


lim Т(а) = 6 Н(2)6 (2-1) gz (H.6) 
a>Q 


QED 

The previous theorem leads to the hypothesis that if 
y is on the closed contour C, anda H(z) is analytic on C and has 
only simple poles inside C, then 


б Н(2) 6(2-р) 92 = 2ujH(v) (н.т) 


where tho integral is taken counter-eloekwise around C. 
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